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Computation  of  Spectral  Data  for  a  Josephson  Junction  Circuit 
E.    G.    Johnson,    Jr.    and  D.    G.    McDonald 

A  computer  program  has  been  developed  to  study  power  flow  between  different 
frequency  channels  in  a  Josephson  junction  circuit.      This  paper  discusses  the  mathema- 
tical assumptions  used  to  get  such  results.      They  are  the  trapezoidal  approximation  from 
spline  theory  and  the  use  of  a  finite  range  of  frequencies  to  characterize  the  frequency 
spectrum.      This  paper  describes  the  program  and  provides  the  FORTRAN  listing,    flow 
charts,    and  discusses  how  to  use  the  program.      A  discussion  of  possible  sources  of 
errors  is  also  included. 

Key  Words:    Differential  equation;  fast  Fourier  transform;  Josephson  junction;  nonlinear 
integral- differential  equation;  spline  theory. 

1.     INTRODUCTION 

The  central  problem  of  interest  is  studying  the  power  flow  between  various  frequency  channels 
in  a  Josephson  junction  circuit  that  is  appropriate  for  harmonic  mixing.     To  be  realistic,    the  circuit 
must  include  shunt  capacitance  and  shunt  resistance  as  well  as  source  andload  impedances  all  of  whichmay 
be  frequency   dependent.     The  general  circuit  we  have  chosen  to  investigate  is  shown  in  figure  1. 

If  the  impedance  Z.,  of  figure  1,  are  not  frequency  dependent,    it  is  relatively  easy  to  solve  the 
differential  equation  for  the  circuit  and  this  solution  is  included  as  an  appendix.      The  result  is  what  we 
term  the  transient  program,    i.  e.  ,    it  gives  the  behavior  of  the  circuit  as  a  function  of  time  starting 
from  specified  initial  conditions.     After  a  few  oscillations,    the  transient  behavior  reaches  a  good 
approximation  to  a  steady  state  oscillation  and  thus  provides  an  estimate  of  a  consistent  set 
of  currents,    of   currents,    voltages,    and  phases    which  we   use   in  the   main  program,    the    steady  state 
program. 

The  steady  state  program  imposes  strict  periodicity  on  the  solution  as  required  for  Fourier 
harmonic  analysis.     Steady  state  solutions  are  generated  in  two  steps,    first  with  the  Z.  frequency 
independent  and  second,    using  the  full  integro-differential  equation  for  general  Z..     In  the  notation  of 


Josephson  junction  when  we  are  given  the  applied  signal  V  .     The  V     is  specialized  to  be 


figure  1,    we  want  the  steady  state  power  for  each  frequency  present  in  the  detector,    Z_,    and  in  the 

.re  given  the  applied  signal  V  .     The  V     is  speci 

V      b    +  bT    cos   (NT  s  +  9T  )  +  bTT  cos   (Ns  +  9TT) 

O  I        O  L  Li  L  U  U  U    J 


V     =  Vjb_  +  bT    cos   (NT  s  +  9T  )  +  bTT  cos   (NTTs  +  9TT)  |   .  (1) 


Here  V     is  the  voltage  scale  factor,   b     is  the  dimensionless  d.  c.    voltage,   b      is  the  lower  applied  r.  f. 

O  O  L 

amplitude,    8L  is  its  phase,    and  b^    9y  are  the  same  quantities  for  the  upper  applied  r.  f.  .     The  integers, 

"frequencies"  N     and  N     generate  multiples  of  the  phase  s,    o<  s  -  W  t  ^  2tt  where  —  is  the  base  period 

Jj  u  o  uu0  c 

of  the  system.     Since  our  strategy  is  to  perform  a  discrete  Fourier  analysis,    obviously  we  must  adjust 

the  discrete  spectrum  to  be  consistent  with  the  applied  frequencies  NT    and  NTT,    the  rotation  frequency 

2eV 
Nj  of  the  self  oscillation  W  =  — —  =  N  U)     and  all  possible  linear  combinations  of  these.     The  b     is 

n  i    o  o 

adjusted  to  allow  only  periodic  solutions  with  each  frequency  related  to  this  base  frequency,  U)  ,    by  an 
integer.     Note  that  bQ  is  a  function  of  b^    h^,   9^,    e        and  N  ,   as  well  as  the  impedances   Z.. 
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This  report  is  subdivided  into  six  additional  sections.     (2)    The  scale  adjustment  is  defined  and 

the  basic  equation  system  is  determined.      (3)  The  approximations  necessary  to  solve  the  basic  equation 

system  are  given.      (4)  The  FORTRAN  Program  and  flow  chart  used  to  solve  the  equation  system  are 

described.    (5)  We  present  a  sample  collection  of  output  from  the  computer  showing  the  following  for  each  case: 

(a)  the  time  plot  of  i   ,   (b)  the  magnitude  of  the  frequency  spectra  for  i   ,    (c)the  magnitude  of  the  frequency 

spectraforV       ,    (d)  the  relative  phase  of  (i    -  V       )  at  each  frequency,   (e)  the  amplitude  spectrum  of -—,  (f)  the 
out  i.        out  as 

steady  state  power  into  the  detector  Z   ,  (g)  the  steady  state  power  across  the  junction  device,   and 
finally  (h)  a  summary  chart  showing  such  things  as  the  scales  of  amperes,  volts,  resistance,  frequency, 
the  b   ,   b    ,   b    ,    9    ,    8       and  the  power  and  other  items  at  the  following  selected  frequencies,   namely 
d.  c.  ,    N    ,    N    ,    the  direct  beat  frequency  N    ,    and  the  rotation  frequency  N  .      (6)  The  summary  gives 
general  conclusions  and  recommendations.      (7)  Finally,    we  have  the  appendix  describing  the  transient 
program. 

2  .     THE  BASIC  EQUATION  SYSTEM 

2 .  1     Definitions 

The  following  definitions  are  noted: 

Z.  are  the  frequency  dependent  impedances  for  the  circuit. 

R  is  the  shunt  resistance  of  the  Josephson  Junction. 

Z     represents  any  residual  impedance  of  the  Junction  device. 

C  is  the  capacitance  of  the  Josephson  Junction. 

cp  is  the  Josephson  phase  function. 

J  is  the  Josephson  critical  current. 

V     is  the  voltage  scale, 
o  & 

W<#,(uu   )'  is  the  radian  frequency  scale. 

R      is  the  resistance  scale, 
o 

A     is  the  current  scale, 
o 

P     is  the  power  scale. 

K  is  the  impedance  operator. 

K      is  the  d.  c.    part  of  the  impedance  operator. 

D  is  the  integral  equation  operator  in  the  circuit  system. 

U  is  the  net  d.  c.    current  needed  to  drive  the  system  at  the  steady  state 

rotation  frequency  N  . 
a    is  the  characteristic  radian  frequency  squared  for  critical  current  in 

dimensionless  form. 

ANU  (v)  is  the  characteristic  frequency  for  loss  due  to  RC  type  effects. 

ANUP  (v    )  is  the  characteristic  frequency  for  loss  due  to  the  rest  of  the 
P 
circuit. 

ANUT  (v _)  is  the  total  loss  characteristic  frequency. 

V,  VA     (v,  v   )  are  the  various  voltages  in  dimensionless  form. 


*(Note  that  the  first  symbol  is  the  FORTRAN  one  and  the  second  is  the  algebraic  form.     This  is  only 
shown  for  the  less  obvious  cases.  ) 


PHA  (cp    )  is  the  phase  with  the  rotation  and  directly  induced  r.  f.    signal 

removed. 
PHB  (cp    )    is  the  directly  induced  r.  f.    signal  phase  plus  the  rotation  and 

the  initial  phase  at  s  =   0. 

PHC  (cp    )  is  the  initial  phase  at  s  =  0. 
c 

THL,  THU    (9    ,  9     )  is  the  phase  of  the  lower,   upper  r.  f.    frequency. 


2.  2.  Equation  System 
We  define 

S  =  iu  t,    V     =  -r-^  ,  A     =  V   /R   ,    Z.  =  z.R   ,    J  =  A   a, 
oo2e  oooiio  o 

C  =  c,/R   U)  ,    R  =  R    r,,  i,  =  A  I,    (<t  =   1  to  6),    and  UO     =  2tt/  T 
loo                 o    1       t,  o  <C  o 


(2) 


Note  that  a,    R    ,    and  (JJ     are  fixed  by  requiring  c    to  be  unity  in  the  dimensionless  system,   by  requiring 

ar     =  —  ,   and  by  using  a     near  1.     This  allows  a  good  initial  guess  for  starting  slope  and  d.  c.    current 

1  1 

in  the  program.     Near  the  hysteresis    region     a1  would  have  to  be  carefully  adjusted  to  allow  study  of 

this  region. 

Further  definitions  are: 


V     =  ■=—  -rr-  -  y    -r-  =  V  V,    i    =  A  a  sm  cp  , 
J       2e    dt  o  ds  o         4  o  T 


dV,  dV,  ...  VT        A       , 

2       „    „        2  .      dV  ,  .  L  o   dcp 

i_  =  C  -r—-  =  V   C  -= —  Ul     =  A    -—  ,    and  l     =  -—  =  -?- 

5  dt  o       ds        o  o  ds  6        R        r      ds 


T       .  V        dVl 

Therefore  l     =  1 .  +  i_  +  l,  =  A       a  sm  cp  +  —  +  -— 

145        6  ol  r  ds  J 

Now  V    s  y  V.,   i  =  i.  +  i_,   i,  =  A  I,,   i.  =  A  I.  , 

s  o    1  1         2       2  o  2       1  o  1 


i  =  A   I   ,    V     ,  =  V   V0  . 
o  o        out  o    2 

Therefore  we  have  the  dimensionless  operator  equations: 


V*3 


'[¥■-]• 


dV 
I,   =  a  sin  cp  +  VI  r.  +  -r—  , 
1  1       as 


I0=(l+z3-1Z2)I1+z3-1V, 


V1  =  Z1I0+  V2    '    and 


V2  =   Z3J2  =   Vl  +  V  • 


(3) 


(4) 


(5) 


(6) 


(7) 


These  are  reduced  to  a  single  operator  form. 

Vj  =  Zj  (1  +  z3"    z2)  I1  +  ZjZ3"    V  +  z2I    +  V,    or 

vi  =  [zi  +  zz+  Vs'S]  [asincp+drJ 

+  {[Z1+   Z2+ZlZ3"lz2]/ri+ZlZ3"1+   !JV- 

We  finally  define 

K=  Zj+  z2+  Sj^"1^.    V*  1/rj,    DHK'1(z1z3"1+    1)  -  ^-  (*f(z°)~  *  +   l)>    and  Vp  ^  (z^0)"1  +   l) 

o  o 

<8) 

as  the  d.  c.    part  of  operator,    where  z.    is  the  D.  C.    value  of  z.,    and 


vT=  v+  vp 


to  get  the  following  pair  of  operator  equations 

-1  dV 

K        (Vn)  =  a  sin  cp  +  — -+  v^V  +  DV,    and 
1  as  1 

v.*.  <9) 

ds 
As  the  reader  can  see,    if  the  z.  are  frequency  independent  then  the  D  operator  is  zero  and  we  have  a 


pure  differential  equation.     If,    as  would  usually  be  the  case,   there  is  frequency  dependence,    then  D 

4  < 

s.,    namely  z     ~  10    ,    z     ~  1,    and  z_ 


4  <  „8 

is  not  zero.    However,    if  we  put  in  estimates  for  the  values  of  z.,   namely  z     ~  10    ,    z     ~  1,    and  z_  ~  10  , 


we  have  D  ~  10      .     If  a  ~  1,   then  our  problem  is  dominated  by  the  differential  equation  and  the  integral 
equation  part  can  be  used  in  an  iteration  procedure.     The  acid  test  is  to  run  the  program. 

As  we  are  interested  only  in  strictly  periodic  solutions,  we  assume  the  following  forms  for  the  Cp 

dcp         dcPb 

and  V.     First  we  split  cp  =  cp    +  cp ,;  V  =  V    +  V,    and  note  that  cp     =Ns+cp    +cp,  and  V,    =  N,  +  -rp-d  s  . 

Ta         b  a  b  bled  b  1       ds  ds 

We  require  cp     to  be  defined  by 

i  d  ^H  dCPH  dCp^ 

K     (viA.c»^  +  \-dT+D^'  (10) 

d   s 

and  cp     to  be  defined  by 
a 

dV 

u  =  — ;-2-+  a  sin  (cp     +  cp,  )  +  v    V    +  DV    .  (11) 

ds  a         b  T    a  a 

Here  u  =  K         b     -  V^N,   . 

o  o  11 

The  purpose  of  this  split  is  to  remove  a  ramp  term  and  to  get  the  first  order  induced  phase  function,    so 

that  the  remaining  part,   cp   V     are  periodic  and  have  a  chance  for  rapid  convergence  to  the  correct  solution. 

Since  V  ,    cp   ,   cp     are  all  periodic  function  of  s,   the  exact  solution  would  be  periodic.     We  use  the 
a        a       d 

Fourier  space  form  to  get  an  idea  of  the  structure  of  D.      Thus  with  -a>  s  t  ^  m, 


cp    (s)  =   ,S  e'lU  cp    (*.),    V    (s)  =  ,2  e"1,t/S  V    {I),    (note  that  V    (0)  =0.) 
a-o  aa-o  a  *  a' 

-iN    s  iN    s  "iNTis  iNTTS 

cpd(s)-e        L5d(NL)+e  9*  (N^)  +  e  V/^+e      U    ^(Njj)  -,  (12) 

.,  2lT  , 

and  D(s)  =  E  e*   S  D(-t).     ("We  have  made  DV  =    f     -j-  D(s  -  s')  V(s').) 

The  equation  (10)  defines  cp     as 

1  bL     _i9L 

K      (NL)_Te  =     ("iNL)+  VT+  °(NL)     ("1NL)  {pdtNL)"  (13) 

The  cp    (N    )  is  gotten  by  U  replacing  L  in  the  above  expression,     cp      is  the  complex  conjugate.     We  note 
d      U 

that  N     and  N      cannot  be  zero. 

U  i-i 

The  equation  (11)  can  be  written  as 


K  ?  ^  d        if 

u6     =     (-it)    +  v(-il)+D(l)  (-U)    cp    M)+af    -^  e1*8  sin(cp    +  cp   )  .  (14) 

o  T  a  I 

•'o 


2tt  """a  "   ^b' 


where  6     is  the  Kronecker  delta.     We  further  note  that  D(-t)  is  given  as 


5(1)  =  (z1(l)z3~1  (l)  +   l)/KM)  -  (Zl(o)z3_1(o)+   l)/K(o), 
where 

KM)  =  z    (t)+  z    (l)+  z    (t)  z  At)  z(l)~l  and 

K(o)  =  K     . 
o 

It    appears  to  be  best  to  do  the  calculation  in  the  time  domain  and  then  do  the  spectral  transform. 
Thus  we  have 


r~      -iNLs 

cpd(s)  =  2  Real      ^(N^)  e 


"iN    s] 


cp,  (s)  =  cp    (s)+  N  s  +  cp   , 
b  d  1  c 


d  /"2TT  d    ' 

u  =  —  V,  +  a  sin  (cp     +  cp    (s))  +  V_V    +/         —  D(s-s')  V    (s') , 
dsa  ah  T    a     •'  o        2rr  a 


(15) 


and 


va  =  a^-  ^a 


as  our  basic  equation  system.     To  progress  further,    we  must  make  numerical  approximations. 


3.     THE  APPROXIMATIONS  NEEDED  FOR  ACCURATE  NUMERICAL  RESULTS 

Wehave  to  choose  a  discrete  set  of  s-points,     N  ,   because  we  must  take  a  finite  amount  of  com- 
putation time  and  must  keep  the  entire  data  set  in  the  machine  to  allow  fast  fourier  transformation. 

N 
The  set  of  points  has  to  be  equidistant  and  equal  to  2    °,   where  N     is  some  integer.     Practical  con- 
siderations restrict  the  range  of  N    from  8  to  11. 

3.  1     Constraints 

A  discrete  set  of  points     implies  frequency  folding,   hence  we  require  the  following  conditions 

to  minimize  this  effect.  .. 

N2 

(a)  The  highest  frequency  of  interest  must  be  ^  —  .     All  frequencies  above  this  value  are  to 

8 

be  treated  as  unreliable,   due  to  the  method  of  analysis.    Computational  noise  must  also  be  small. 

(b)  Our  demand  accuracy  for  analysis  is  near  5%.  We  make  checks  on  the  results  by  increas- 
ing N  as  the  occasion  demands.  If  the  spectrum  does  not  change  to  the  demand  accuracy,  our  analy- 
sis is  treated  as  adequate  and  accurate  to  that  order.     Incidently,    if  the  value  of  U  is  correct  to  5% 

or  less,   then  the  lower  frequency  spectra  is  also  correct  to  that  accuracy. 

(c)  This  5%  computational  accuracy  level  suggests  using  N_  =  256  or  512  with  an  occasional 

extension  to  2048  points  for  special  studies.  ,. 

2 

(d)  We  allow  z.(-t-)  to  have  arbitrary  values  for  0  s  l^— —and  require  the   remaining  t  values  to  be 

fixed  by  z.*M)  =  z.(N     -  I). 

(e)  In  D(-t)  and  cp  ,(<£,),  we  use  inplace  of  (-i-C)  the  periodic  expression(e  -  l)/h,  where  h=  2rr/N     is 

d  <2 

the  interval  between  points  in  the  time  space.     This  allows  folding  and  preservation  of  reality  condi- 
tions.    The  sensitivity  to  this  change  is  checked  by  the  doubling  of  points. 

(f)  Because  D(s  -  s1)  is  a  small  correction  to  the  differential  equation  system,   we  replace, 

with  no  qualms,   the  integral  by  a  sum  using  the  trapezoidal  rule.      Thus 

N2-l 

rds'D(s-s')V(s')=     E     rr  D(j  -  j')  V  (j1)-  (16) 

/    ~ST  a  j'=0N2  a 

•'o 

dV 

(g)  To  complete  the  approximations,   we  must  decide  how  to  approximate  — —  type  terms. 

Finite  differencing  does  not  allow  high  enough  precision  in  the  fit.     Since  the  spectral  information  for 

dV  2 

power  uses  — —  in  its  current,   we  need  to  have  smoothness  in  the  results  to  at  least  order  h   . 

, 

Additionally,   we  want  to  avoid  absolute  instabilities  and  permit  the  application  of  quasilinear  theory 

to  get  the  b     value  needed  to  allow  strict  periodicity.     These  requirements  make  it  necessary  to  use 

°  2 

the  spline  method     using  at  least  the  trapezoidal  rule.     We  assume  the  trapezoidal  rule  everywhere  in 

the  calculation.     The  results  of  this  analysis  should  be  smooth  functions  (accuracy  order,   h    ),    hence 

a  frequency  separation  of  true  results  from  noise  due  to  the  discrete  analysis  of  the  system  and  an 

2 

inherent  overall  accuracy  of  h  .     Periodicity  should  tend  to  force  uniform  error  throughout  the  time 

interval.     Sharp  spikes  in  the  time  dependence  of  cp     wm  negate  the  uniformity  of  the  error  and  make 
higher  frequencies  spectral  information  less  accurate. 

(h)    Now  the  problem  is  split  into  a  part  (a)  which  neglects  D  and  a  part  (B)  that  takes  D  into 
account. 


3.  2    Part  (a) 


Here  we  have  the  simplified  system  of  equations  using  the  trapezoidal  rule. 


cp   (j+  D  =  cp   (j)  +  T     V   (j+  1)+  V   (j) 


va(j+  D  =  va(j)+  -  |q,(j+  D+  q 


|[va(j+i)+va(j)] 

q^)  =  -a  sinlcpa(j)+  cpb(j)      -  VTVa(j)+  u. 


V    (o)  =  V    (N,),    and  cp    (o)  =  cp    (N). 
a  a      2  a  a      2 


The  u  is  adjusted  to  allow  this  result.     The  eq  (19)  can  be  rewritten  with  definitions 


hj  =  (2  -  hvT)/  (2  +  hvT)  ,   h2  =  -  ah/  (2  +  hvT)  ,   h3  =  2h/  (2  +  hvT), 


The  boundary  conditions  are 


Va0  +   1)  =  h1  Va(j)+  h3u+  h     I  sincp(j  +  1)+  sincp(j) 


(17) 
(18) 
(19) 


We  generate  two  solutions  starting  with  cp    (o)  fixed  and  unchanged  for  a  given  calculation  run.     To  avoid 

numerical  noise  problems,    experience  shows  that  it  is  absolutely  necessary  to  fix  the  cp    (o),   when  u  is 

a 

adjustable. 


3.  3    Calculation  Sequence 


<P    (j  +   1 )  -  cp    (j) 

a  »      a 


of  the  n+  1  iteration  changes  by 


(1)  We  iterate  for  a  given  j,   j+  1  pair  until 
less  than  EPI  =10". 

(2)  We  compute  two  quasilinear  solutions  using  the  following  starting  conditions 

6cp  (1)  (o)  =  0  ,    6V    (1)(o)  =  1  ,    and  6cp    (2)(o)  =  6V    (2)(o)  =  0. 
a  a  a  a 

In  case  one  6u  =  0  and  in  case  two  6u  =   1.      The  two  equation  sets  are,    with  q    (j)  =  cos   cp(j), 


(i)  ( 

6cp    l  '  (j  +  1)  =  Sep    v 
a  a 


and 


(i) 


(i), 


6V  0+   D=  h    6V    x  '(j)+  h.  6 

a  la  5     c. 


i,(J)  +  |[5Va(1,(j+l)+6Va(»U)], 


q2(j+   D6cp(i)  0+   D+  q2(j)  6cp( 


%]  ■ 


(20) 


These  can  be  rewritten  as 


6cp  (j  +   1)  =  6cp    (l)" 

a  a 


i)(j)  +  |[6va(i)a+i)+6va(i)a)]. 


6Va(i)(j+   l)=Tt(i)/T2(i) 


with 


?l(i)  *    [hl  +  h2lq20+   UJ  6Va(l)^)+  h3621+  h25cP(l)<J>  [q20+   D+  q20)l.  (2D 


(i)  h 

and  T^'-       l-h2lq20+l)      . 

(3)    We  use  the  values  computed  at  cp    (N   ),    6cp      (N  ),    etc.   to  satisfy  the  following  require- 
ments  and  hence  fix  P  and  6u: 

cp   (N  )  -  cp   (o)+  0   6cp    (1)  (N,)+  6U  6cp    (2)(N,  )  =  0, 
a      2  a  a  2  a  2 

and  V   (N   )  -  V    (o)+  0(6V    (1)(N   )  -   1)+  6u  6V    (2)(N   )  =  0.  (22) 

a      2  a  a  i.  a  2 

The  new   cp    (j)  ,    V    (j)  »   an<*  u  are  now  defined  by 
a  a 

cp    0)=  cp    0)   (old)  +  0  6cp    (1)(j)+  6u6cp    (2)(j)  , 
a  a  a  a 

V    (j)  =  V    (j)    (old)  +  P  6V    (U0)  +  &u  6V    (2)(j)  .  (23) 

a  a  a  a 

and  u  =  u    (old)    +  6u  . 

In  the  computer  program  the  corrections  to  the  initial  guesses  of    u    and  V    (o)  are  bounded  because 

initially  the  first  guess  may  not  be  in  the  circle  of  convergence  except  when  the  computer  program  in 

the  appendix  is  used.     Thus,    it  is  necessary  to  have  built  into  the  calculational  procedure  a  means  to 

search  for  the  region  where  the  quasilinear  technique  will  work.     The  computer  program  has  allowed 

jumps  in    6u    and  in  6  V    (o)  that  shrink  automatically  as  each  computation  is  made  with  a  new    u    and 
a 

V  (o)  .     If  after  ten  iterations  the  convergence  has  not  taken  place,   then  the  V    (o)   is  set  equal  to 

3.  3. 

V  (N  )  using  the  idea  that  it  is  more  likely  to  be  compatible   with  a  periodic  solution  because  the 
transient  should  be  less.     Then  we  repeat  the  above  quasilinear  process  and  the  window  shrinking 
process   for  up  to   eight  times.     If  the  convergence  is  particularly  difficult,   the  shift  in  V    (o)  can 
take  place  up  to  ten  times.     Finally,   a  relaxed  test  for  acceptance  of  the  results  as  converged  is  made 
so  that  a  spectral  analysis  can  be  made. 

3.  4    Problems 

There  are  two  main  sources  for  convergence  problems:     (1)  The  number  of  computed  points 

is  too  small  when  there  are  sharp  spikes  in  the    V      and  -: —    .     This  causes  extensive  computational 

noise.      (2)  The  initial    V    (o)  and  u  are  badly  guessed.     Both  problems  are  particularly  bad  when 

there  are  strong  B    ,    B      terms  and/ or  when  a,   v      are  very  large  numbers. 
V         J-i  T 


3.5    Part  (P) 

3 
Here  we  use  the  full  quasilinear  theory.        First,    we  define 

N2-l 

q    (j)=  -a  sincp(j)  -      £     D(j  -  j')  s    V    (j ')  +  u  -  vT  V    (j) 

J     =  0  2 

q2(j)  =  coscp  (j) 

q30)  =  h3H   (Va^+   1)_  Va0)>-  °' 5  (^ i  <J  +   D+  q^))'  (24) 

q40)=  hji  (cpa(j  +    1)  -q>a(j»  "  0-5  (Va(j+   D+  Va(j))| 

We  note  that   q_,     q      are  the  errors  due  to  not  fitting  the  equation  system  to  a  quadratic    spline.     We 

use  these  terms  to  define  the  error  signal  in  the  quasilinear  theory.     In  addition  to  the  terms  needed 

to  correct  the  error  signal,    we  must  add  effects  due  to  needed  changes  in  V    (o),   and  u  to  allow  pe- 

a 

riodic  solutions.     Since  the  quasilinear  theory  integrates  over  the    entire  domain,   the   u  is  sensitive 

to  the  global  solution,    and  periodicity  tends  to  force  global  accuracy;  we  should  retain  a  global  ac- 

2 
curacy  of  order  h   ,    and  a  smooth  spectrum  to  that  same  order.     This  should  be  accurate  enough  to 

separate  the  true  solution  from  the  effects  of  the  finite  analysis.     Again,    spiking  will  make  the  higher 
frequencies  less  accurate. 

The  quasilinear  differential  equations  are  replaced  by  the  trapezoidal  rule.     The  basic  equa- 
tions become 


.       ( 
6cp 
a 


and 


°G  +   D=  6cpa(l)0)+|["6Va(i)(j)+  6Va(i)(j+   1)1  -  q4(j+   D  *\ 
6V    (i)(j+   D  =  h.  6V    (l)(j)+h ,  [q,(j)6^(l)(j)+  q,(j+  D  6cP    ^O  +   1)1  (25> 

+  h36x2  -  q3(j+   1)6^    ,    for  i=  1,2,3. 

The  initial  conditions  are  6cp    ^'(o)  =  0,  for  all  i;  the  6v    l1'  =  0  for  i  =  1,  2,   and  the  6v         (o)  =  1. 

a  a  a 

The  periodic  conditions  at  j  =  N     requires 

6cp(1)(N   )  +  6u  6cp(2)(N   )  +  3  6cp(3)(N2)  =  0, 

6 V(1)(N  )  +  6u  6V(2)(N  )  +  P  (6V(3)(N  )  -   1  )  =  0.  (26) 

c.  2  « 

We  use  these  equations  to  fix  6u  and   P.     We  compute   for  the  different  a,   a  new  cp     = 

old  cp    +  a(6       cp    +  6u6        cp    +  P6       cp    )   ,    and  a  similarily  constructed  new  V     which  are   then  inserted 

a  a  aa  aN2-!22 

to  compute  a  new   q      and   q      set  for  each  a.     The  sum  of  these;    namely,    P(a)  =      N^     (q       +  q      ), 

3  4  jVo 

gives  us  an  error  signal.     Our  choice  for  the  correct  cp    ,    V      is  fixed  by  the  requirement  the  correct 

a         a 

a   has;  P(a)  <P(o).     The  smallest  P(cp)  improves  the  calculation  the  most  and  guarantees  convergence. 
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The  a  chosen,    is  then  used  to  construct  the  new  cp     and  V     and  the  entire  part  (3)  quasilinear  analysis 

a  a  2 

is  repeated.     This  process  is  continued  until  10  iterations  have  been  performed  or  until  P(o)<  2h   . 

The  results  of  this  process  constitutes  the  best  accuracy  that  can  be  obtained  for  a  global  solution 

subject  to  the  inherent  errors  of  the  trapezoidal  method. 

Next  we  compute  by  fast  Fourier  transform  the  desired  spectra.     This  completes  the  discussion 

of  the  assumptions  made  to  get  numerical  results.     The  next  section  has  the  FORTRAN  program  listing 

and  a  flow  chart  of  the  basic  program  flow.     In  particular,   we  show  the  input  data  requests. 

4.  THE  FORTRAN  PROGRAM  AND  FLOW  CHART 

Here  we  have  listed  the  Fortran  programs  that  give  the  main  details  of  the  computation.     The 
microfilm  and  fast  Fourier  transform  programs  are  not  listed.     The  former  depends  on  the  computer 
installation  and  the  latter  should  be  available  and  optimized  at  each  computer  installation.     After  the 
Fortran  list,   figure  2,    a  flow  chart,    figure  3,    shows  the  bulk  data  and  computational  flows  between 
programs. 

5.  THE  SAMPLE  COLLECTION  OF  OUTPUT 

This  section  discussed  the  printed  and  the  microfilm  output  for  a  single  computed  case.     The 
purpose  is  to  give  general  information  to  allow  use  of  the  Fortran  program  by  an  interested  party.     No 
attempt  is  made  to  report  all  the  possible  data  that  exists  as  a  consequence  of  the  computations. 
Selected  information  generated  as  a  consequence  of  this  program  will  be  described  in  separate  reports. 
To  use  the  program  correctly  in  detail,    one  must  read  the  Fortran  listing. 

5.  1     The  Printed  Output 

Figure  4  shows  the  printed  output  from  a  single  computed  case.     The  first  row  of  data  gives 
the  number  of  data  points  used,   the  number  (11)  corresponds  to  2048  points  in  time  and  1024  different 
frequencies.     The  J  (1.  -4)  is  the  supercurrent  in  amperes.     The  Al  (1.  25-0)  gives  the  scale  adjustment 
defined  on  page  4.     C  (8.  2277-15)  is  the  capacitance  in  farads.     R  (1.   ±  1)  is  the  number  of  ohms. 
EB1    (1.  -8)  is  the  demand  accuracy  for  each  step  interation     and  scales  the  overall  demand  accuracies 
in  thevarious  iterations  of  part  (a)  and  part  (3).     The  EP6  (1.  -0)  can  be  set  to  increase  demand  accuracy 
of   part  (3).     EP6  equal  1.    uses  the  values  stated  on  page  10  for  the  converged  exit.     EP6  less  than 
one  increases  the  accuracy  demand  and  EP6  greater  than  one  decreases  accuracy  demand. 

The  first  Zl,    etc.  ,    row  gives,   with  19  equal  1,   the  complex  values  of  Zl,    Z2,   and  Z3  for  all 
impedances  as  if  they  were  frequency  independent.     Subsequent  rows  of  Zl  type  data  show  specific 
frequency  exceptions  to  the  frequency  independence  assumption.     These  are  overlaid  on  the  previous 
values.     Please  note  that  19  is  one  unit  larger  than  the  corresponding  frequency  number,   thus   1  is  the 
d.  c.    case.     In  our  sample  output  we  have  changed  the  first  harmonic  of  the  Z3  impedance  from  (1.  +  8,  0.  ) 
to  (1.  +  5,  0.  ).     The  others  have  remained  unchanged. 

We  continue  the  row  by  row  analysis  of  the  printed  output.     After  the  impedance  spectrum  is 
defined,    we  next  specify  the  lower  frequency,    NL(4),   the  upper  frequency,    NU(5),   at  which  the  two  r.  f . 
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PROGRAM  JOSPHC 
C      THIS  PROGRAM  COMPUTES  FREQUENCY  SPECTRA  FOR  A  CIRCUIT  MODELING 
C      THE  POINT  CONTACT  JOSEPHSON  JUNCTION.  THE  ANALYSIS  USES  A  TIME 
C      PERIOD  WITH  N2  POINTS.  NO  IS  THE  ORDER  2**N0=N2   NO.LE.  11  . 

DATA  IP 1=3.141 592653 59) .( IM=(0.,1. ) ) , ( PHS=2.06778E-15 > 

DIMENSION  ZK1025)  »Z2(1025)  »Z3(1025) 

COMMON  /5/  BI»S1»S2»S3 
1»A1 

COMMON  /16/  V16.P16.H8 

TYPE  COMPLEX  Zl ,Z2 »Z3 .K ,DT ,ZF »ZH»ZG »KL »AN2» ATMP.KU.ATN.AN6.AN7.ZI 
1»CBU»CBL  »ZD 

TYPE  REAL  KO 

COMMON  /DAT/  PI.IM.PHS.U 

DIMENSION  K(1025> ,DT(1025> .KL(3),KU(3> 
1   »I2T(1025> »VAT<1025> 

EQUIVALENCE  (K(  1) ,DPH(1 ) ) » ( DT ( 1 ) ,DVA(1>) 

COMMON  /l/  PHC»VA(2050).DPH(4»2050>» 

1DVA( 4.20  50) »H,H1»H2.H3.H4,EP1.EP2.EP3.LBX(5> ,LBTR(5) .EP4(6) 

EQUIVALENCE  (DPH(6151)  .ZKlD  *  (  DVA  (  6151 )  *Z2(  1 )  )  »  ( Q4<  1 )  »Z3(  1 )  ) 

COMMON  /2/D(2050)  ,PHA(2050)  .PHB(2050>  »QK2050)  *Q2(2050)  .03(2050)  » 
1Q4(20  50>  ,KEY(6) .NO ,N1 ,N2 »NH ,NP ,P0 ,V0» AO ,V»   THL .THUtBO.BL »BU.NL ,NU 
2.NB.KL.KU.  AN5.H5,   ANU , ANuP > ANUT »A .Rl .RO ,C. WO.T ,   R, PI2 »AN2 ,AN4» 
3AN3»MM.N5»N6.N7»ISKPD»SC»  OUT < 6. 11 ) »NR.    ICONV.AJ. ICASE. 

4POW1.POW2.POWJ.POW I N.K0.EP5.H7 
5.CBU.CBL 

EQU I VALENCE <  X(l) ,DPH(1 ) )»(Y(1) »DPH(2051) ) . ( VAT( 1 ) »DVA( 1 ) ) » 
1(  I1T(1  )  .DVA (2051)  )  .(VOUT(D  .DVA (4 101)  )  »(X1(D  ,DPH(  1025)  ) 

TYPE  COMPLEX  VAT. I  IT »VOUT. I2T  .IM 

DIMENSION   XQ025)  .  Y  ( 1025  )  »  I  IT  (  1025  )  »VOUT(1025)  ,X1  (  1025  )  »Y1  (  1025  ) 

EQUI VALENCE (Y1(1),DPH( 4101) ) , ( I2T ( 1 ) .Q3( 1 ) ) 

1  READ  2.N0.AJ.A1.C.R.EP1.EP6 

2  FORMAT* I5.6E12. 4) 

PRINT  200.N0.AJ.A1.C.R.EP1.EP6 

200  FORMAT(*  NO » J .Al »C »R »EPl .EP6  */I5»6E12.4> 

N2=  2**N0  tNH=N2/2   $AN3=N2   SPI2=  2.*PI  $AN4=PI2/AN3  $N5=N2+1 

N6=N2+2   $N7=N2-1   $AN2=-IM*AN4  $MM=N0~1  $SC=.5/AN3  $NP=NH+1 

XT=PI2*C   S  R7=PHS/(XT*AJ)   $  A=R7/ ( R*A1*. 1 )**2 

R0=SQRTF(A*R7)  $T=R0*XT  *W0=PI2/T    $NR=NH/2 

Rl=  R/RO   J  ANU=1./R1   $VO=PHS/T  $A0»V0/R0  $PO=VO*AO   $AN8=Nh 

EP2=EP1*A   $EP5=.1*A 

S1«R1*SQRTF(A)  $S2=A1/10.   $  S3=EXPF ( .2-Sl ) 
C      SET  I8  =  I9«=1  FIRST  TIME  TO  GENERATE  ALL  THE  SAME  Z1.Z2.Z3. 
C      TO  ALLOW  SELECTED  FREQUENCY  DIFFERENCE  USE  18=0  19  THE  PARTICULAR 
C      HARMONIC.  TO  EXIT  USE  I8=-l  AND  19  THE  LAST  PARTICULAR  FREQUENCY 
C      HERE     I9=NP  AND  Z1.Z2.Z3  ARE  REAL  AS  MUST  BE  Zlll)  ETC. 
C      19=1  IS  D.C.  CASE  ADD  1  TO  HARMONIC  ORDER  TO  GET  CORRECT 
C      LOCATION  IN  DATA  STORE. 

5  READ  3.I8.I9.ZK 19) .Z2( 19) »Z3( 19) SPRINT199. 19 »Z1 ( 19) .Z2 ( 19) »Z3(I9) 

3  FORMATf I3»I5»3C(E12.3»E12.3) ) 
IF( 18)4,5.6 

6  DO  7  I=2»NP  $ZKI)=ZKlJ  $Z2(  I>=»Z2(1>  $Z3(D=Z3(1) 

7  CONTINUE      $G0  TO  5 

4  DO  8  1  =  1, NP  $Z1( I)=Z1(  IJ/RO  $Z2 ( I ) =Z2 ( I ) /RO  $Z3 ( I ) =Z3 ( I ) /RO 
ATMP  =  1.+Z1(  D/Z3( I)  SK( I )=Z1( I >+Z2< I )*ATMP  $DT (  I  )  =  ATMP/K( I ) 

8  CONTINUE    SANUP=DT(1)  $ANUT=ANUP+ANU   $D0  9  I=2.NP 

9  DTI  D=DT(  I)-DT(l)  $KO  =  K(D    $  DT(1>=0. 
READ  10. NL.NU.ISKPD. ICASE 

C      NL  IS  LOWER  APPLIED  FREQUENCY   NU  IS  UPPER 

10  F0RMATUI5) 

C      ISKPD=+1  DO  THE  D  STRUCTURES.  ISKPD=0  SKIP  D  STRUCTURES  HERE  D=0. 
PRINT  201.NL. NU.ISKPD. ICASE 

201  FORMATJ*  NL »NU, ISKPD. ICASE  *»4I5> 

C      ICASE  IS  INCREMENTED  TO  LABEL  EACH    COMPUTED  CASE  LOAD  LAST  VALUE 

TF(NL-NU) 11.12.12 
12     PRINT  202»NL.NU      SCALL  EXIT 

202  FORMAT (  *  NL  TOO  LARGE  *  .215) 
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11     IF(NU-NR>13»14.14 

14  PRINT  203.NU.NR 

203  FORMAT (  *  NU  TOO  LARGE  *   215) 
13     NB=NU-(NU/NL)*NL 

KEY(1)=0   SKEY(2»=NB  $KEY(3»=NL  $KEY(4)=NU   $KEY<5>=NH 
IF( IS<PD>15»19 

15  DT(NP)=.5*DT(NP)  $ZH=CEXP ( AN2 )  $ZG=CONJG(ZH>$DO  16  J-1.N2  $D(J)*0« 
ZG=ZG*ZH  $ZF=1.  $DO  17  I=2»NP  $ZF=ZF*ZG  $ZD=DT<I)*ZF 

17  D( J)=D(J)+REAL(ZO) 

16  D(J)=D(J»/AN8      $  DT(NP)»2.*DT(NP) 
READ  18.IDRED  $60  TO  ( 19.20 >♦ IDRED 

C      DECIDE  PRINT  OF  D  IDRED=2  YES   -1  NO 

20  PRINT  204»<D(J) ,J»1.N2) 

204  FORMAT  (  *  D  SET  *  /(10E12.3>> 

19     CONTINUE         $H=AN4   $EP3=2.*H*H*EP6 

H4=H*.5   $  HT=H*ANUT  S  HU=1./  ( 2.+HT  >  $H3=2.*H*HU  $H2«=-A*HU*H 
H1=(2.-HT)*HU    »H9=H4*H2  $  H7=HU 

AN7=  AN2*flU    $ATMP=(CEXP<AN7)-1.)/H  $ATN=( ATMP+ANUT+DT(NU+1 ) ) 
KU=  0.5/(ATMP*K(NU+D*ATN)  $KU(  2  >=ATMP*KU  $KU(  3  )  =ATMP*KU<  2  ) 
AN6=AN2*NL   $ATMP= ( CEXP ( AN6 > -1. ) /H  $ATN= ( ATMP+ANUT+OT (NL+1 ) ) 
KL  =  .5/(ATMP*K(NL+1)*ATN)    $KL  (  2  )  =KL*ATMP   SKL  (3  )=ATMP*KL<  2  > 

C      Nl  IS  RAMP  FREQUENCY  RELATIVE  TO  BASE   USE  PHC=-PI  IS  Nl  NOT  ZERO 

C      PH0=PHC  AT  S=0.   BL.THL  IS  LOWER  APPLIED  RELATIVE  VOLT  AND  PHASE 

C      BU  »THU  ARE  THE  UPPER 
CALL  CHARTA 

C      NOTE  BL=K0*BL  USED  IN  TRANSIENT  PROGRAM.  SAME  FOR  BU.  HERE  Ko  IS 

C      Zl/RO.   HERE  R0=R*A1*A/10. 

2022  READ  2 1 *N1 .PHC.BL, THL»BU»THU 
BL«  KO*BL   $  BU=KO*BU 

C      ADD  EACH  TIME  ICOuNT  »1  IGUS=-1  FOR  SMALL  POWER  INCREMENTS. 
C      ADD  ICOUNT  =1  IGUS=1  AND  USE  RESULTS  OF  TRANSIENT  PROGRAM. 
C      ADD  ICOUNT  =1   IGUS=0   TO  START  NON  RF  GUESS  SEQUENCE. 
C      ADO  ICOUNT  =0.-1  TO  GET  NEW  POWER  VALUES 
2020   READ  2041 » ICOUNT» IGUS.BI »V16»P16 
2041   FORMAT(2I5.3E20.9) 
IF(EOF»60>22»23 

22  CALL  EXIT 

23  IF( ICOUNT12022. 2022,2023 

2023  AN5=H*N1   $  KEY(6)=N1 
S4=N1*S2 

PRINT    2111»N1»PHC,BL»THL»BU»THU 
2111       FORMAT(*    Nl.PHC. BL.THL. 8U.THU    *    I5»5E12.4) 

21  FORMAT!  I5.E20. 9    /2E20.9/2E20.9     ) 
CBU=BU*CEXP(-IM*THU)       $    CBL=BL*CEXP <-IM*THL > 
ZF=CEXP(AN6)     SZG=CEXP(AN7)    $ZH=CONJG< ZF)*KL*CBL 
ZI=CONJG(ZG)*CBU*KU  *DO    24    I=1»N5    $AI=I~1 
ZH=ZH*ZF    $ZI=ZI*ZG    $    E*2»*(ZH+ZI) 

24  PHBU)=E+PHC+AI*AN5  $ICONV  =  30  $ICASE=  ICASE+1 
IF( IGUS>302»300»303 

302  BI=U    $V16=VA(1)     $P16=PHA(D     $    GO    TO    303 

300  BI=A*S3*(SQRTF(1.+S4*S4>-S4> 

301  V16=P16=0. 

303  CONTINUE 

PRINT  2113  .ICOUNT. IGUS.BI, V16.P16 
2113   FORMAT(*  ICOUNT . IGUS »BI .V16 .P16  */2 15 .3E20.9 ) 
CALL  PARTA    $IF  (ICONV)  25.26.26 

25  PRINT  205*  Nl.ICASE 

205  FORMAT (  *  NO  CONVERGENCE  FOR  Nl=  *  15.  *  AND  CASE  NUMBER-  *  1 5  ) 
GO  TO  2020 

26  ICONV=10   SCALL  PART3   $IF( ICONV > 27 .31 .31 

27  PRINT  205.N1.ICASE 
31     CALL  CHARTB 

GO  TO  2020 

18  FORMAT (15) 

199    FORMAT!*  Zl .Z2.Z3 ( 19 )=  *  15 »3C< E12.3 .E12.3 > > 
END 
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SUBROUTINE  SCALE ( 19. 199 > 
C      PURPOSE  IS  TO  CONVERT  THE  DC  DOWN  BY  10-7  AND  FORM  LOGlO(Yi) 
C       THE  AMPLITUDE  OR  POWER  SPECTRA  5  DECADES  ONLY. 

TYPE  REAL   MIN.MAX 

COMMON  /4/  ICLE(8.2»6>  »IDC(8) 

COMMON  /l/  E(6151) » Yl < 1025 > »F< 11299 > 

GO  TO  (1.2.3.4.5.6)199 

3  IDC(3>=ICLE( 19+1.1. 3>=ICLE( 19+1. 2. 3 > =0   $  RETURN 

1  CONTINUE 

2  CONTINUE 

4  CONTINUE 

5  CONTINUE 

6  CONTINUE 

C      POWER  CASE  AND  VOLT  CASE 

7  MIN=0.  $  MAXaO. 
IF(  19)8,9 

9      IDC(I99»=-7   $  Y1(3)=1.E-7*Y1(3) 

8  DO  11  I»1.1024 
IF(Y1( I>)12»11»13 

12  IF(MIN-Y1(I> >11.11.14 

14  MIN  =  YK  I)  $  GO  TO  11 

13  IF(MAX-Y1(I> )15, 11.11 

15  MAX  =  Y1U) 
11     CONTINUE 

IF(MIN)31.32.32 

31  IMIN=AL0G10<-MIN)  i    MI N=-l.* < 10. ) **IMIN  SICLE (  19+1 .1 » 199 )  =  IMIN-5 

32  IF(MAX»33»33.34 

34     IMAX=AL0G10<MAX)$MAX=U0.)**IMAX  SICLE ( 19+1 .2 » 199 )= IMAX  -5 

33  DO  21  1=1.1024   $  I F( Yl ( I ) ) 22 »21 .23 

22  Yl( I)=-(ALOG10(Y1( I )/MIN)+5. )   $  IF ( Yl < I » ) 21 .21 .25 
25     Yl( I)=0.  $  GO  TO  21 

23  Yl( I)=ALOG10(YK I)/MAX)+5.   $  IF( Yl ( I )) 24.21.21 

24  Yl( I »=0.   $  GO  TO  21 

21     CONTINUE   $  RETURN  $END 


SUBROUTINE  PART 

DATA(  ALPHA=0..1...2».1..01..005) 

DIMENSION       ALPHA (6> 

DIMENSION  21(1025) .Z2 ( 1025 ) »Z3 ( 1025 > 

COMMON  /5/  BI, SI .52. S3 
1.A1 

COMMON  /16/  V16.P16   .H8 

TYPE  COMPLEX  Zl .Z2 .Z3 .K ,DT .ZF.ZH.ZG.KL .AN2.ATMP.KU.ATN.AN6.AN7 ,ZI 
l.CBU.CBL 

TYPE  REAL  KO 

COMMON  /DAT/  PI.IM.PHS.U 

DIMENSION  K<1025>  »DT(1025) .KL(3) .KU(3> 
1   .I2TU025)  »VAT(1025) 

EQUIVALENCE  (K(l) ,DPH(1) ) . ( DT ( 1 ) ,DVA ( 1 ) ) 

COMMON  /l/  PHC,VA(2050> .DPH( 4.2050 ) , 

1DVA(4.2050) .H.H1»H2.H3»H4.EP1»EP2.EP3.LBX(5J .LBTR ( 5 ) .EP4 ( 6 ) 

EQUIVALENCE  (DPHC615D  »ZK1>)  »  (  DVA(  6151 )  »Z2(  1  )  )  .(Q4(l)  .Z3«D  ) 

COMMON  /2/D(2050)  »PHA(2050)  .PHB<2050>  .QK205  0)  »Q2(  2050  )  »Q3  <  2050  )  . 
1Q4(2050),KEY(6) .NO .Nl »N2 .NH.NP.PO.VO.AO.V »   THL .THU.BO »BL »BU,NL »NU 
2.NB.KL.KU.  AN5.H5,   ANU.ANUP .ANUT .A »R1 *R0 .C.WO, T .   R.PI2 »AN2 ,AN4» 
3AN3.MM.N5.N6.N7.ISKPD.se.  OUT ( 6» 11 ) »NR .    ICONV.AJ. ICASEt 

4P0W1.P0W2.P0WJ.P0W I N.K0.EP5.H7 
5.CBU.CBL 

EQUIVALENCE!  X(l) ,DPH( 1 ) ) . < Y( 1 ) »DPH(2051) ) » < VAT ( 1 ) .DVA( 1 ) ) ♦ 
K I1T(1) .DVA(2051) ) » ( VOUTI 1 > .DVAI4101 ) ) »<X1(1) »DPH(1025)) 

TYPE  COMPLEX  VAT . I IT .VOUT. I2T  »IM 

DIMENSION   X  ( 1025  ).Y(  1025)  »  I  IT  (1025)  .VOUT  (1025)  .XI  (1025)  .YK1025) 

EQU I  VALENCE  (YK1)  ,DPH(4l0l)  >  .  (  12  T(  1  »  »Q3(  1  >  ) 

DIMENSION  Tl<4)  .T2<4)»T3«4> 
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ENTRY  PARTA 
C      THIS  PART  USES  PURE  DIFFERENTIAL  EQUATION  AND  TRAPAZOIDAL  SPLINE 
C      APPROACH  TO  GET  PHA.VA  SHAPE  IN  TIME 

DO  1  I-1»N6    $VA< I )=0. 

1  PHA(I)=0. 
EP77=EP5 

EP78=.01*A  $ZP=EP78*.3 
EP74=EP78 

IF(N1>2.3 

3  U=  A*SINF(PHB<1>)  SGOTO  483 

2  U«BI 
VA<N5)=V16 

483    ICONVl=12 

DO  460   119=1,10 

BET«6AM«XI»0. 

PHA(D=P16  $  VA(1)=VA<N5)   *  VA(N5)=0. 

ICONVl=ICONVl-2   $  IF( ICONV1 ) 61 »61 .620 
61     ICONV=2  $  GO  TO  4 
620    ICONV=ICONVl 

4  ICONV=ICONV-l  $IF( ICONV ) 460 .461 »461 

460  CONTINUE   $  IF< ABSF ( PHA ( N5 > -PHA ( 1 > )-100.*EP2 >411 ,413 »413 

411  IF(ABSF(VA(N5)-VA(1) )-100.*EP2 > 412 .413 »413 

412  ICONV=l  $  GO  TO  405 

413  ICONV=-l$  GO  TO  405 

461  U=U+GAM 
VA<1)=VA( 1)+XI 

TA=H3*U   $TB  =  PHA(D+PHB(1)   $Q3  ( 1  '  =SINF(  TB) 
401    Q2<l>=COSF(TB)  $DPH(  1  »  1  )=DPH(  2  »1  )=DVA(  2.1  )  =  0« 

DVA(1.1>=1.   $  T3<2)=H3  ST3<1>=0. 

DO    6    1  =  1. N2       $13=1  +  1       $TD=H1»VA(I  >+TA-trH2*Q3II)         $    IV=5 
8  TB=PHA( I3)+PHB( 13)     SQ2 ( I  3 > =COSF I TB)    $    TC=1 .-H5*Q2( I3> 

IV=IV-1 

IF(    mio.ii.ii 

10  PRINT  13    SCALL  EXIT 

13  FORMAT <*  CONVERGENCE  PROBLEMS  IN  PART  A*  ) 

11  Q3( I3»=SINF(TB) 

VA( I3>=TD+H2*Q3( I3» 

TB=PHA(I)-PHA(I3)+M4*(VA(I3)+VA( I> ) 
DEL=TB/TC 

PHA( I3>=PHA( I3)+DEL         $        DEL=ABSF(DEL ) 
IF(DEL-EP1)7»8»8 
7      TB=PHA( I3>+PHB( 13) 

Q3(I3>=SINF(TB>$TC=    Q2  ( I  >  +Q2  (  13  )  $VA<  13 )=TD+H2*Q3 <  13) 

Q2( I3>=COSF<TB> 

TB=H5*Q2(I3»         $DO  5  11=1.2 

Tl(  I1»=DVA(  I1.I)*(H1+TB>+T3(  I1»+H2*DPH(  I1,I>*TC 

T2(I1)=1.-TB 

DVA(  I1.I3)=T1(ID/T2<I1) 

5  DPH(I1.I3»  =  DPH( I1,I)+H4*(DVA< 1 1 » 13 >+DVA( II » I ) ) 

6  C9NTINUE     $GO  TO  12 

12  B1=PHA(1)-PHA(N5)  $B2=VA< 1 > -VA< N5 ) 

1450   A11=DPH(1,N5)  $A12=DPH( 2.N5 ) 

A21=DVA(1»N5)-1.  $A22=DVA<2»N5> 

DET=A11*A22-A12*A21  $B1«B1/DET  $  B2=B2/DET 

XI  =A22*Bl-A12*B2   $GAM=A11*B2-A21*B1 

IF (ABSF( GAM  J -EP5» 17,17.19 
17  IF(ABSF(XI)-EP74)14,14.19 
19     GAM=SIGNF(EP77.PHA( 1)-PHA(N5» ) 

IF(A8SF(PHA(1)-PHA(N5> ) -.1 ) 199. 199.198 

198  XI=SIGNF(ZP»VA(4>*(VA(1)+VA(4)-VA(3)-VA(2) »> 
GO  TO  200 

199  XI=SIGNF(EP78.(VA(1)+VA(4)-VA(3)-VA(2) )*VA(4) ) 

EP78=.9*EP78 

200  EP77=.90*EP77 

14  PRINT  1000»U.VA(1).PHA(1).GAM.BET.XI 
1  »VA(N5>.PHA(N5» 

1000   FORMAT!*  U.VA.PHA  *  /(6E13.4>> 

I F ( ABSF ( GAM  > -EP2 ) 400 .4.4 
400    IF(ABSF«XF)-EP1)40.4»4 
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40  RETURN 

405    PRINT  999»<PHA< 17) ,VA( 17) ,I7=1,N5> 
999    FORMAT(8E12«4) 

ENTRY  PARTB 
C      THIS  PART  GETS  THE  INTEGRAL  EQUATION  PART  USING  QUASILINEAR 
C      ANALYSIS  AND  THE  TRAPAZOIDAL  METHOD  FOR  THE  ERROR  ESTIMATE 

DPH(1»1)»DPH<2»1)=DPH<3,1>=DVA(2»1  >=DVA«  1*1  )=0.  $DVA<3*1)=1. 
60     J6*l  SGO  TO  44 

30     GO  TO  (21*28,28, 28, 28, 24*24, 60»35>IWAY 

21  A12=-H4$A11=A22  =  DET1=  1.  $  T2 ( 2  )  =T2 ( 3  )  =0.  $T3(2)=H3   $T3(3)»0. 
DVA(2»1>=DPH<2»1>=0. 

DO  22  1=1. N2    $12=1+1  $A21=-H2*Q2( 12)  $DET=DET1-A12*A21 

T3(1>=-Q3(I2) 

T2(l>=-Q4< 12)  $DO   23  16  =  1*3 

B1=(DPH(  16*  I)+H4*DVA(I6*D+T2(  16  >  )/DET 

B2=(H1*DVA( 1 6. 1 )+T3( I6>+H2*Q2( I )*DPH( 16*1 ) J/DET 

DPH(I6»I2>=   B1*A22  -B2*A12 

DVA(I6»I2>=   B2*A11  -B1*A21 

23  CONTINUE 

22  CONTINUE 

Bl=  -DPH(1»N5>   $B2=  -DVA(1»N5) 

All=DPH(3»N5>  $A12=DPH(2»N5> 

A21=DVA(3*N5)-1.  $A22=DVA( 2 »N5 > 

DET=A11*A22-A12*A21  $B1=B1/DET  $  B2=B2/DET 
XI  =A22*B1-A12*B2   $GAM  =  AH*B2-A21*B1 
199=2   $  GO  TO  42 

24  TEMP=EP4( 1)SI99=1$IC0NV=IC0NV-1  $J6=1$IWAY=8     SIF ( ICONV) 39,26,26 
26     DO  250  L6=2»6   $  IF ( EP4( L6)-TEMP ) 25 »250,250 

25  TEMP=EP4(L6>  $J6=L6 

250    CONTINUE    $  IF  (J6-D  62.62*28 

28  DO  32  I=1»N2 

DPH(2»I»=    ALPHA (J6>*  DPH ( 1  *  I >+PHA ( I > 
32     DVA(2»D=    ALPHA(  J6>*  DVAd.I  )+VA(I)     $TC=U+ALPHA(  J6)  *GAM 
GO  TO  <46»49> 199 

29  EP4(J6)=0.  $D0  34  I=1»N2 

SMALL=Q4( I)**2+Q3( I )**2   $  IF(SMALL-EPl  )341*341.344 

344  PRINT  345»SMALL»U*I,J6 

345  FORMAT(*  ERROR  AT  POINT  *  »2E12. 3*215) 
341    EP4< J6)=EP4( J6)+SMALL 

34  CONTINUE 

TEST=EP4(J6>   $  IWAY=J6   $  J6=J6+1 
IF(  TEST    -EP3)51»30.30 

51  IWAY=9  SGO  TO  46 

62     PRINT  63»EP4»EP3    SRETURN 
39     PRINT  41,EP4,EP3 

41  FORMAT(  *   10  LOOPS  IN  QUASI  NO  CONVERGENCE  */(6E13.4)) 

35  RETURN 

42  DO  43  I=1*N2 

DPH(1»I)=DPH(1,I)+XI  *DPH(3»I )+GAM*DPH(2» I > 

43  DVA(1»D=DVA(1,I>+XI  *DVA  (  3  »  I  )+GAM*DVA(  2  »  I  > 
IWAY=J6=2  $  GO  TO  30 

44  DO  45  I=1*N2 
DPH(2»I,=PHA(  I)  $DVA<2»D=VA(  I) 

45  CONTINUE    $  TC=U  SGO  To  49 

46  DO  47  I=1»N2   $PHA< I )=DPH( 2 » I )  $VA( I ) =DVA ( 2 ♦ I ) 

47  CONTINUE       $U=TC  $GO    TO    30 

49  DO       52     1=1. N2       $TB=DPH( 2 » I> +PHB ( I )     SQ2 ( I > =COSF( TB) 

Ql(  D=-A*SINF(TB)+TC         -ANUT*DVA(  2  .  I ) 
IF( ISKPD)54»52 

54  DO      53    I3=1»N2 
L=I-I3 
IF(L)55.56.56 

55  L=N2+L 

56  CONTINUE 

53  Ql(  D=Q1(  I)-D(L  +  1)*DVA(2»I3) 

52  CONTINUE 
DVA(2»N5)=DVA(2.1) 
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DVA(2»N6)=DVA<2»2) 

DPH(2*N5)»DPH(2»1J 

DPH(2.N6)=DPH(2»2) 
63  FORMAT (*    NO    CONVERGENCE    POSSIBLE      *    /(6E13.4)) 

Q2<N5)=Q2<D       SQ1(N5)=QKD    $Q1 (N6 ) =Q1 (2 ) 

DO    50    I=1»N2      $    12=1+1    $13=12+1 

04 ( I2)=DPH(2.I2>-DPH(2»I>-H4*(DVA<2»I2>+DVA(2»I> ) 
50  03(  I2)=H7*(DVA(2»I2)-OVA(2»n-H4*(QK  I2>+Q1(  I)  )     ) 

Q3<1>  =Q3<N5)$Q4(1>  =Q4(N5»  SGO  TO  29    SEND 


SUBROUTINE  CHART 

DIMENSION  Zl ( 1025 ) »Z2 « 1025 )»Z3( 1025) 

COMMON  /5/  BI.S1.S2.S3 
1.A1 

COMMON  /4/  ICLE(8»2»6) »IDC(8) 

COMMON  /DDC/LU»LUC*IFL 

COMMON  /DD/  IN.IOR»lT,IS»IC»ICC»IX.IY 

COMMON  /TAB/  VI ( 60 ) , V2 < 121 ) »V3 ( 19 ) 

TYPE  COMPLEX  Zl .Z2 »Z3 .K »DT .ZF .ZH.ZG.KL »AN2 .ATMP.KU.ATN.AN6. AN7.ZI 
l.CBU.CBL 

TYPE  REAL  KO 

COMMON  /DAT/  PI.IM.PHS.U 

DIMENSION  K < 1025 ).DT (1025) »KL(3)»KU(3) 
1   »I2T(1025> .VAT(1025) 

EQUIVALENCE  <K( 1 ) »DPH( 1 > > » (DT < 1 ) ,DVA( 1 > ) 

COMMON  /l/  PHC,VA(2050» »DPH( 4*2050 ) . 

1DVA(4»2050)»H,H1,H2,H3*H4,EP1,EP2.EP3*LBX(5),LBTR(5),EP4(6) 

EQUIVALENCE  (DPH(6151) *Z1(D ) » ( DVA< 6151 ) »Z2 ( 1 ) ) , < Q4( 1 ) ,Z3( 1 > > 

COMMON  /2/D( 2050 >»PHA( 2050' »PHB(205Q) »Q1 < 2050 > ,Q2( 2050 ) ,Q3(2o50)* 
104(2050) »KEY(6) ,N0 *N1 *N2 »NH»NP ,P0»V0»A0,V,   THL,THU»BO,BL*BU*NL,NU 
2»NB»KL»KU.  AN5.H5,   ANU,ANUP»ANUT »A»R1 »R0 ,C» W0,T ,   R.PI2 ,AN2 ,AN4, 
3AN3»MM»N5»N6»N7»ISKPD,SC.  OUT ( 6» 11 ) »NR»    ICONV. AJ, ICASE. 

4P0W 1 » POW2  » POW J » POW I N  » KO , EP5 . H7 
5.CBU.CBL 

EQU I  VALENCE  (  X  ( 1  >  »DPH<  1 )  )»(Y(D  .DPH(2051>  )  »(VAT(D  ,DVA(1»)  . 
1(  I1T(1)  »DVA(205D  )  »(VOUTU>  »DVA(4101)  )  »  (XI  ( 1 )  ,DPH(  1025  )  ) 

TYPE  COMPLEX  VAT* I IT »VOUT* I2T  »IM 

DIMENSION   X(1025) ,Y<1025) » I  IT (1025) .VOUT(1025) ,X 1(1025) »Y1(1025) 

EQU  I  VALENCE  (YK1»,DPH(  4101)  )  »  (  I2T  (  1 )  »Q3(  1 )  ) 
ENTRY  CHARTA 

REWIND  1 

WRITE  TAPE  1*Z1,Z2»Z3 
REWIND  1 

IN=IOR=IT=IC=0  $IS=1  $LU=LUC=0  $IFL=30 
CALL  DDINIT(3»20H  EGJOHNSON.EXT  3706  ) 
S1»S1/(R1*ANUT> 

PRINT  100»A1»S1   $  PRINT  110»A»ANUT 
110    FORMAT  (*  A»ANUT  *,2'E12.4) 

ENCODE ( 40* 100* VI )A1»S1 
100    F0RMAT(16H  Al »SQR( A) /NUT • /  »2E12«4> 

IX=89  $  IY=500  $  CALL  DDTAB  $  CALL  DDTEXT15.V1) 

ENCODE(480.1»Vl)NQ,N2»NB»NL»NU.ANU»ANUT,A.V0.A0,PQ.WQ»RO»C»R,H,AJ» 
IT 
1      FORMAT (89H  THE  JOSEPHSON  SPECTRUM  CONTROL  PARAMETERS  THAT  REMAIN 

1C0NSTANT  '/FOR  THE  ENTIRE  RUN.*/   .40H  NO, DATA  POINTS.BEAT ,LOWER»U 
2PPER  FREQ.'/  »5I5»  48H'/  NU.NU  TOTAL.SCALED  A.VOLTS. AMPERES*  WATTS 
3'/    »  6E12.4*  62H  '/RADIAN  FREQ, OHMS*  CAPACITANCE.  RESISTANCE*  ST 
4EP  INTERVAL'/   .  5E12.4.36H  '/SUPER  CURRENT  AMP. .PERIOD  SEC.'/ 
5   2E12.4.24H  '/CASES  TO  FOLLOW.'.  > 

ENC0DE(32»99»LBX(1) ) 
99     FORMAT (32H  JOSEPHSON  CURRENT  IN  TIME  CASE       ) 

IX=89$IY=950  JCALL  DDTAB  $CALL  DDTEXT ( 60.V1 )  SCALL  DDFR  SRETURN 
ENTRY  CHARTB 
READ  TAPE  1»Z1.Z2»Z3 
REWIND  1 
C      Q2(D  IS   A*SINF(PH> 


Figure  2.   (Continued) 
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DO  2  I=1.N2 

2  Q2t I»«A*SINF(PHA( I)+PHB( I) ) 
X1(D=0.  $DO  3  1=2.256 

3  Xl(  D»XK  I-D+4.  $IFL=60 
ENCODE* 8 .98. LBX (5) ) I CASE 

98     FORMAT* 15) 

NS  =  N2/257+l  500  4  I8=1»NS  $I9  =  256*( 1 8-1 >  $00  5  1  =  1*256  $11=1  +  19 

5  Y1(I)=Q2(U»   $CALL  FILMGRAF  (  XI  » Yl  *256.LBX.  1*0  ) 

4  CONTINUE 

DO  6  1=1. N2 

6  Q2<  I>=Q2<  D+QK  I>+ANU*VA«I> 
DO  7  I=1»NH   $X(I)=Q2« I ) 

7  CONTINUE  $D0  8  I=NP»N2 

8  Y( I-NH)=Q2( I)    $CALL  REORDER ( X.Y.MM )  $CALL  CFFTRC(X. Y.MM.SC. 1 ) 
CALL  REALTRAN(X.Y.MM.l.l)  $D0  901  =  1. NP  $  I  IT (  I ) =CMPLX(X< I ) »Y< I ) ) 

90     CONTINUE 

C         CALL  REORDER(A.B.M) 

C         REVERSIBLE  PERMUTATION  OF  REAL  SEQUENCE 

C         FROM  FIRST-LAST-NORMAL  SEQUENCE 

C         TO  ODD-EVEN-REVERSE  BINARY  SEQUENCE. 

C         SEQUENCE  LENGTH  IS  N  =  2**M 

C         CALL  REALTRAN(A.B.M.NE.INV) 

C         IF(INV.GT.O)  UNSCRAMBLE  THE  TRANSFORM  OF  A  REAL  SEQU'-NCE. 

C         IF(INV.LT.O)  SCRAMBLE  THE  TRANSFORM  OF  A  REAL  SEQUENCE. 

C         INPUT  AND  OUTPUT  ARE  IN  NORMAL  SEQUENCE, 

C         SEQUENCE  LENGTH  IS  N  =  2**M 

C         NE  MUST  AGREE  WITH  SIGN  OF  EXPONENT  IN  TRANSFORM  DEFINITION. 

C         INNER  LOOP  SINES  AND  COSINES  COMPUTED 

C  RECURSIVELY  BY  SINGLETONS  2ND  DIFFERENCE  ALGORITHM, 

C  INITIALIZED  FROM  A  DATA  TABLE. 

C         DISCRETE  COMPLEX  FAST  FOURIER  TRANSFORM. 

C         CALL  CFFTRC(A.B.M.SC.NX) 

C         INPUT  A(J)  +  I*B(J)  IN  REVERSE  BINARY  SEQUENCE. 

C         OUTPUT  A<K>  +  I*B(K)  IN  NORMAL  SEQUENCE. 

C         SEQUENCE  LENGTH  IS  N  =  2**M 

C         SC  IS  REAL  SCALING  MULTIPLIER. 

C         NX  IS  THE  SIGN  OF  THE  EXPONENT  IN  THE  TRANSFORM  DEFINITION. 

C         INNER  LOOP  SINES  AND  COSINES  COMPUTED 

C  RECURSIVELY  BY  SINGLETONS  2ND  DIFFERENCE  ALGORITHM. 

C  INITIALIZED  FROM  A  DATA  TABLE. 

I1T(1>=  UT(1)+N1*ANU   $IlT(NU+l)  =  IlT(NU+l,  (  KU(  2  )  *ANU+K.U<  3  >  )  *CBU 

1 1 T ( NL+1 > = I  IT ( NL+1 )  +  < KL ( 2  > *ANU+KL ( 3  > > *CBL 

DO  10  1=1. NH  $X(I)=VA( I ) 

10  CONTINUE  $D0  11  I=NP.N2 

11  Y(I-NH)=VA( I)  $CALL  REORDER ( X.Y.MM )  $CALL  CrFTRC(X. Y.MM.SC . 1 ) 
CALL  REALTRAN(X.Y.MM.l.l)  $DO  12  1=1. NP  $VAT ( I )=CMPLX(X ( I ) . Y( I ) ) 

12  CONTINUE 

VAT(1>=VAT(1)+N1   $VAT(NU+1>=VAT(NU+1>+KU(2>  *CBU 
VAT(NL+1)=VAT(NL+1)+KL(2)   *CBL 

VJUN=V0*VAT(1)   $  AIJUN=A0*I1T(D  $  ENCODE*  64  .9111,  VI  >VJUN,AUUN 
9111   FORMAT ( 30H  V JUN ( VOLTS ) .  IJUN< AMPERES M /    »  2El7«5> 
DO  13  1=1. NP 
VOUT(I)=  Z2< I )*I1T( I )+VAT( I)    $I2T ( I ) =VOUT( I ) /Z3 < I ) 

13  CONTINUE   $D0  14  1=1.1024     $X1(I)=I-1 

14  Y1(I)=0. 

C      DROP  NP  POINT  IN  THESE  PLOTS.  SEE  PRINT  OUT  FOR  IT. 
POW1=POW2=POWJ=0. 

NS  =  NH/129  +  1  $D0  15  18=1, NS   $  I  9=128* ( 18-1 )   $K7=I8-1 
DO  21  1=1,128   $  I2=(I-l>*7+3   $11=1  +19 

21  Y1(I2)=CABS(I2T(I1) )   $ENC0DE(40,97,LBTR( 1 ) > ICASE, 18 
97     FORMAT(26H  AMPLITUDE  I2TRAN.  ICASE=  215.  4HPART   ) 

CALL  SCALE(K7.1) 

CALL  FILMGRAF(X1.Y1.1024»LBTR.1.0» 
C      PLOT  OF  CURRENT  .12.  AMPLITUDE  THROUGH  DETECTOR  Z3.  AS  A  FUNCTION 
C      FREQUENCY 

DO  22  1=1.128  $I2=( I-l>*7+3  $11=1+19 

22  Yl( I2>=CABS«V0UT( ID )  $ENCODE( 40.96.LBTR* 1 ) > ICASE. 18 
96     FORMAT (26H  AMPLITUDE  VOuTTR.  ICASE=    215.  4HPART  ) 
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CALL  SCALE! K7. 2) 

CALL  FILMGRAF(X1»Y1»1024»LBTR.1»0) 
C      PLOT  OF  OUTPUT  VOLTA6E  ACROSS  DETECTOR 
DO  23  1=1.128  $I2=(I-l)*7+3   $11=1+19 

23  Y1(I2>=  CANG( I2T( ID )-CANG(VOUT( II) > 
ENCODE(40»85.LBTR(l)  )  I  CASE.  18 

85     FORMAT (26H  PHASE  I2-VOUT  TRAN  ICASE=  215.  4HPART   ) 

CALL  SCALE(K7»3> 

CALL  FILMGRAF(X1»Y1»1024.LBTR.1.0) 
C      PHASE  OF  I2-VOUT 

DO  24  1  =  1.128  $I2=(I-D*7+3  $11=1  +  19 

24  Y1(I2)=  CABS(VAT( II) )  $ENCODE (40,95 ,LBTR( 1 >) ICASE. 18 
95     FORMAT (26H  AMPLITUDE  VPHTR.  ICASE=    2l5»  4HPART   ) 

CALL  SCALE(K7.4) 

CALL  FILMGRAF(X1.Y1»1024»LBTR.1»0' 
C      PLOT  OF  AMPLITUDE  OF  VOLTAGE  ACROSS  JUNCTION. PH( DOT) 
DO  25  1=1.128  $I2=(I-l>*7+3  $11=1+19 
TEP=  2.*REAL<  I2T( 1 1 )*CONJG( VOUT ( II > ) )  $POW2=POW2+TEP 

25  Y1(I2)=TEP   $  IF(I9)26»27 

27  TEP=  Y1(3)  =  .5*YK3»   $POW2=POW2~TEP 

26  ENCODE(40»94.LBTR(1> )ICASE»I8 

94     FORMAT (26H  POWER  I2V0UT  SSTR.  ICASE-  2I5»  4HPART     ) 

CALL  SCALEIK7.5) 

CALL  FILMGRAF(X1»Y1»1024.LBTR.1.0) 
C      DETECTOR  POWER  FOR  EACH  FREQUENCY.  P0W2  IS  TOTAL  POWER. 

DO  28  1=1.128  $I2=( 1-11*7+3  $1 1=1+19 

TEP=2.*REAL( I  IT ( II )*CONJG ( VOUT ( II ) ))  $POW1=POW1+TEP 

28  Y1(I2>=TEP  $IF(I9)29.30 

30     TEP  =  Y1(3»=.5*YK3)  $P0W1=P0W1-TEP 

29  ENC0DE(40»93.LBTR<1> )ICASE»I8 

93     FORMAT (26H  POWER  I1VOUT  SSTR.  ICASE=   215.  4HPART    ) 

CALL  SCALEIK7.6) 

CALL  FILMGRAF(X1»Y1»1024»LBTR.1.0> 
C      POWER  IN  22+JUNCTION    POW1  IS  TOTAL  POWER  IN  POWJ  IS  TOTAL  IN  JUN 
C      KO  IS  THE  EFFECTIVE  IMPEDENCE.   NOTE  ALL  PLOTS  ARE   LOG10  EXCEPT 
C      PHASE  CASE.   THE  POSITIVE  .UP  CASE  IS  SCALED  INDEPENDENT  OF 
C      THE  NEGATIVE  CASE, DOWN  CASE.   THIS  IS  DONE  IS  SCALE  SUBROUTINE. 
C      FILMGRAF   THEN  DOES  A  SCALE  CHANGE  TO  MAKE  THE  DATA  FULL  SCALE 
15     CONTINUE 

OUT ( 1 ♦ 1 )  =  ( U+N1*ANUT ) *K0  $OUT ( 2  » 1 ) =OUT (5.1) =OUT ( 6. 1 ) =0. 

0UT(3»D=BL  $0UT(4.D=BU  $OUT  (  3»2  )  =THL  $OUT(  1  »2  )  =PHC+PHA(  1 ) 

0UT(4»2'=THU  $OUT(2»2)=OUT(5»2)=OUT(6»2)=0. 

0UT(1.3)=0UT(1,D*(  UT(1)  +  I2T(D)   $OUT(  2  .3  )  =OUT  (  5  .3  )  =OUT(  6.3  )  =0. 

OUT (3.3) =BL*REAL ( ( 1 1 T ( NL+1 >  + 1 2  T ( NL+1 > ) *CEXP ( I M*THL ) ) 

0UT(4.3>=BU*REAL( ( I 1T( NU+1 >+I2T (NU+1 ) ) *CEXP( IM*THU) ) 

POW I N  =OUT ( 1 »  3 ) +OUT ( 3  » 3  > +OUT (4.3) 

P0W1=P0W1+REAL(  I1T(NP)*C0NJG(V0UT(NP) ) ) 

P0W2=POW2+REAL(  I2T ( NP) *CONJG (VOUT (NP) ) ) 

POWJ  =  REAL(  I1T(1)*C0NJG(VAT(1)  ) +1  IT (NP)*CONJG( VAT(NP) ) ) 

DO  32  1=2. NH 
32     P0WJ=P0WJ+2»*REAL( I  IT ( I ) *CONJG( VAT ( I ) > ) 

DO  34  14=1.6   $  I5=KEY( I4>+1 

OUT( 14,4) =  2.*REAL(I1T( 15 )*CONJG( VOUT( 15) ) ) 

OUT(  14, 5)=  2.*REAL( I2T( 15 )*CONJG( VOUT( 15  >  >  > 

0UT(I4.6)=  REAL(Z1(I5)>   $OUT ( 14.7 )=AIMAG(Z1 ( 15 > ) 

0UT(I4.8)=  REAL(Z2(I5))   $OUT ( 14,9 >=AIMAG(Z2( 15 ) ) 

0UT(I4,10)=REAL(Z3( 15) )   $OUT ( 14, 11 ) =AIMAG(Z3(  15 ) ) 

34  CONTINUE 

DO  35  14=1,6,4   $  DO  36  17=4,5 
36     OUT( I4»I7)=.5*0UT(I4.I7) 

35  CONTINUE 

ENCODE (968  ,9,V2 ) POWJ.P0W1 ,P0W2 .POWIN.KO. ( KEY ( I ) » (OUT( I.J) »J=1»11) 
1.1=1.6) 
9      F0RMAT(28H  POWJ .P0W1 .POW2.P0WIN .KO, • /    5E12.3,2H«/     . 
138HKEY.B.PHC.P0WIN.P0W1.P0W2.,/Z1,Z2»Z3»/        . 
2  6(I4.5E12.3.2H«/,6E12.3»2H«/   )) 

ENCODE (152. 91 »V3>EP1.EP2»EP3.EP4 
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91     FORMATUOH  ERROR  LIMITS  EPl»EP2»EP3  •/  EP4(6)  •/    »3E12.3»2H«/ 
16E12.3.  2H«.  ) 

IX=89   $  IY=950   $  CALL  DDTAB  $  CALL  DDTEXK 121»V2 ) 

IX=89  $  IY=600  $  CALL  DOTAB  $  CALL  DDTEXT(8tVl» 

IX=89   $IY=500     $CALL  DOTAB  SCALL  DDTEXT ( 19»V3 ) 

IX=89  $IY=400  $  ENCODE(480»911.V2> (IDC(I) *I=1*6) »( H ICLE(I.JtL) »I= 
1  1.8) »J=1»2> »L=1.6) 

CALL  DDTAB  $  CALL  DDTEXT (60. V2 )  $  CALL  DDFR 
911    FORMAT (58H  SCALE  EXP  ONLY  DC  SET  THEN  MIN.MAX  SET  FOR  EACH  GRAPH 
1'/   »  6I4*2H«/  »6( 16I4»2H»/) ) 

RETURN 

END 
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(  IF(I8)   )4.5.6 
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(IF(ISKPD))15,19 
15(   READ    )IRED 


19(CONTINUE  ) 

(   CALL  )CHARTA 

2022(   READ  )  ,N1 

2020(   READ  )ICOUNT .  IGUS 


CALL  PARTA 
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(  GO  TO  ^)2Q2Q 


._l 


THE  MICROFILM  PLOT 
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FOURIER  TRANSFORM 
TO  GENERATE  THE 
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THE  INITIAL 
GUESS  FOR  PH 
AND  THE  QUASI- 
LINEAR  PROCEDURE 
FOR  THE  SELF- 
CONSISTANT  PH. 

ENTRY  PARTA 
ENTRY  PARTB 


Figure  3.      The  flow  chart  for  Josephson  steady  state  program. 
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signals  are  to  be  applied,   whether  we  want  to  compute  the  integral  equation  part,    (ISKP  is  zero 
means  that  the  frequency  spectra  of  impedance  is  to  be  ignored.     ISKP  equal  1  means  that  it  is  consi- 
dered and  the  integral  equation  is  now  used.);  and  finally  what  the  initial  case  number  label  for  this 
computation  run  is  to  be,   ICASD(5).     The  remaining  cases  are  incremented  by  one  for  each  new  case. 

The  Nl,   PHC,    etc. ,    row  is  repeated  for  each  new  case.     Here  Nl  (4)  is  the  chosen  rotation 
frequency  for  the  Josephson  oscillator,    PHC  (-3.  1416-0)  is  the  chosen  initial  phase  assuming  PHA 
is  zero  initial  value.     BL  is  the  corrected  value  for  the  lower  frequency  r.  f .   applied  voltage  as  defined 
on  page  1.     This  BL  (2.  4+5)  is  changed  from  the  BL  (7.  68+2)put  into  the  program  and  used  in  the  tran- 
sient program  in  the  appendix.     The  lower  r.  f .   phase  is  THL(0.  )  in  this  example.     The  BU(0.  )  is  like- 
wise scaled  from  BU(0.  )  and  finally  we  have  the  upper  initial  phase  THU(0.  ). 

The  next  data  row  has  two  control  variables  ICOUNT  and  IGUS  as  well  as  input  data  that  is  used 
as  initial  guesses  for  U,    V    (1),    and  cp    (1).     The  ICOUNT(l)  has  two  options--if  it  is  zero,    it  causes 

3.  3. 

the  program  to  request  a  new  set  of  Nl,    PHC,    BL,    THL,    BU,    THU;  and  if  it  is  equal  to  one  it  causes 
the  program  to  proceed  to  compute  the  case  just  read  in.     The  IGUS  has  three  options--if  it  is  -1,   it 
causes  the  program  to  use  the  computed  previous  values  of  U,    VA(1),    PHA(l)  as  the  initial  guesses 
for  the  new  set  of  data;   (this  option  is  useful  if  you    are  simply  making  small  changes  in  previous  data 
cases),    if  IGUS  is  zero,    it  uses  some  computed  guesses  based  on  the  approximation  that  the  r.  f.   power 
is  not  importantly  affecting  the  initial  data,    hence  the  d.  c.    case  can  be  a  good  basis  for  initial  data 
guesses;  and  finally,   if  IGUS   is  one,    it  uses  the  currently  read  in  values  BI  for  U,    Vl6  for  VA(1),   and 
Pl6  for  PHA(l).     These  values  are  gotten  from  the  program  in  the  appendix  or  from  other  calculations. 
The  final  rows  of  information  in  the  printed  output  depends  on  how  the  calculation  proceeds.     If 
the  calculation  is  normal,    then  you  get  what  is  shown  in  the  sample  case.     If  the  calculation  does  not 
work  normally,   then  you  will  get  error  messages  and  some  detailed  printouts.     Please  look  at  the 
Fortran  listings  to  understand  these  messages  and  printouts.     The  U,    VA,    PHA  sequence  is  repeated 
until  there  is  convergence  in  the  PARTA  subroutine.     The  data  printed  out  is  the  current  U(-2.6+l), 
VA(1)  (4.96),   PHA(l)  (-3.62,   6U  (1.6177),   change  in  cp   (1),  the  change  in  VA(1),  the  VA(N5)  (4.6127), 
and  finally  PHA(N5)  (-3.  6492).     Please  note  that  convergence  occurs  when  the  VA(1)  equals  VA(N5) 
and  PHA(l)  equals  PHA(N5).     Also  note  that  the  phase  PHA(l)  is  unchanged  for  the  entire  calculation 
sequence.     Here  (1,  N5)  is  equivalent  to  (0,N   )  in  the  analysis. 

5.  2    The  Microfilm  Output 

We  look  at  selected  plots  of  the  microfilm  output  of  a  single  computer  run.     This  data  corres- 
ponds to  the  same  run  as  the  printed  output.     Figure  5  has  information  that  remains  constant  through- 
out the  different  runs.     The  only  information  that  may  not  be  self  evident  is  the  SQR(A)/NUT  term  (.  5) 
which  is  the  characteristic  frequency  ratio  of  the  circuit.     It  isWa/vrr,  and  it  measures  the  capacitance 
level;   the  smaller  it  is  the  less  important  is  the  capacitance  term  in  the  circuit. 

Figure  6  has  8  plots  of  256  points  each  of  a  single  period  of  computation.     Please  note  that  the 
FILMGRAF  routine  will  independently  scale  the  individual  plots  so  that  the  variations  are  full  scale 
in  the  y  axis.     If  the  range  of  y  variables  are  both  negative  and  positive  then  the  scale  is  plus  and 
minus  one.     If  the  y  variables  are  only  one  sign  then  the  range  is  zero  to  minus  one  or  one  to  zero 
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Figure  5.     Sample  microfilm  output  -  the  control  parameters. 
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depending  on  the  sign.  We  emphasize  that  it  is  possible  for  the  different  plots  to  have  very  different 
scales  for  their  respective  data.  In  this  sample  case  the  figure  6  set  has  exactly  the  same  scale  for 
all  8  plots. 

Figures  7  and  13  are  the  first  two  plots  of  a  series  of  8  for  the  amplitude  at  each  frequency  of 
the  current  through  the  detector.     The  first  figure  gives  the  first  128  harmonics,    the  latter  figure  gives 

the  next  128  harmonics.     The  remaining  768  harmonics  are  not  shown.     The  d.  c.   line  amplitude  is 

-7 
reduced  by  10       in  the  first   plot.     Then,    in  this  set  of  128  harmonics  as  well  as  each  succeeding  set 

of  128  harmonics,   the  largest  positive  value  is  found,   the  l°gin  is  taken  and  the  integer  part  is  extracted. 

Then  all  positive  nonzero  y  values  are  scaled  relative  to  this  integer.     Only  the  top  five  decades  of 

log       amplitude  values  are  computed  and  put  into  the  y  axis.     The  remaining  y  values  that  are  less  than 

this  allowed  5  dacades  are  set  equal  to  zero  and  is  thus  plotted  that  way.     For  those  plots  which  have 

negative  nonzero  y  values  the  largest  negative  value  is  extracted,   the  integer  part  is  removed,   these  y 

values  are  scaled  down  by  that  number,    and  the  log       is  taken  and  those  numbers  that  are  5  decades 

less  than  this  negative  maximum  are  set  to  zero.     The  numbers  are  then  made  negative.     Both  the 

negative  and  positive  values  that  have  been  processed  as  above  are  now  subject  to  further  scaling  in 

the  subroutine  FILMGRAF.     That  scaling  is  the  same  as  already  discussed  for  the  time  series  plots. 

The  above  scaling  takes  place  for  all  the  figures  7  through  18  as  well  as  those  plots  that  have 

not  been  shown.     Figure  19  has  at  its  bottom  a  set  of  scaling  information  which  in   conjunction  with  the 

information  given  at  the  bottom  of  each  figure  permits  the  relative  determination  of  the  amplitudes  for 

each  spectral  line  that  has  been  plotted.     The  row  of  -7's  with  a  zero  mean  that  only  the  phase  differ- 

-7 
ence  series  of  plots  does  not  scale  the  d.  c.    frequency  line;  the  rest  are  scaled  by  this   10      .     This  was 

done  to  allow  detailed  plot  information  of  the  higher  frequencies  without  the  d.  c.    values  causing  scaling 
problems.     After  the  -7  row  are  six  rows  of  numbers.     The  first  eight  numbers  of  each  row  defines 
the  scaling  used  on  the  negative  numbers  in  the  spectrum.     Note  that  only  the  last  row  which  corres- 
ponds to  the  power  spectra  across  the  junction  has  negative  values.     The  last  eight  values  correspond 
to  the  scaling  of  positive  numbers.     The  rows  are  related  to  the  plots  as  follows. 

Row  Plot 

1  Amplitude  of  the  current  across  detector  (AMPLITUDE  I2TRAN) 

2  Amplitude  of  the  Voltage  across  detector  (AMPLITUDE  VOUTTR) 

3  Phase  difference  between  detector  current  and  voltage 

(PHASE  I2-VOUT  TRAN) 

4  Voltage  across  R  (AMPLITUDE  VPHTR) 

5  Power  across  detector  (POWER  I2VOUT  SSTR) 

6  Power  across  Josephson  junction  (POWER  I1VOUR  SSTR) 

We  illustrate  the  reading  of  the  plots  to  get  an  absolute  number  for  the  power  across  the 
junction  at  the  4,  8  harmonics  and  at  the  132,    133  harmonics  respectively.     The  first  two  are  gotten 
from  figure  12  and  the  latter  from  figure  18.     The  conversion  to  watts  is  given  in  figure  5  as    (4.  8830  -11). 
The  fourth  harmonic  has  value  full  scale  (5.  918).     The  eighth  harmonic  has  value  (4.  7  X  (4.  4410)/ 6.  2 
=  3.  34). 
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Figure  6.      The  next  8  plots  give  one  complete  period,    each  plot  256 
points  in  time. 
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Figure  6.     (Continued) 


27 


JOttPMSON  CUMCNT   IN  Tl?€  CASC 


8  000*000 
J  Ml. MO 


t  wooes 


-2  K0O02 
•2  940*002 


2  140O92 

2  HO. 002 


»*  !f?2  «•»  OS  ;ltli  «t  !?:M-M 


.l.iK  ».  3 


2K  »••»!•.  r^flM 


Figure  6.    (Continued) 
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JOSEPHSON  CURRENT  IN  TIME  CASE 


N»(  t  iflt 

sc«:# 


Hiit 

0.000*000 
0.000*000 


Ha. 

1.020*003 

1.020*003 


•   Hi* 

-2.560*002 
-2.560*002 


2.560*002 
2.560*002 


B«i>   !5"*2  H«y  03    1124)   at   17:50:53  Ll.L2«   I.    0 


256  ^a<»t».   ►"'••• 


Figure  6.      (Continued) 
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JOSEPHSON  CURRENT  IN  TIME  CASE 
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Figure  6.      (Continued) 
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JOSEPHSON  CURRENT  IN  TIME  CASE 
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Figure  6.      (Continued) 


31 


JOSEPHSON  CURRENT  IN  TIME  CASE 
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Figure  6.      (Continued) 
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JOSEPHSON  CURRENT  IN  TIME  CASE     6 
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Figure  7.     The  first  frequency  block  for  the  LOGjq  amplitude  of  the 

current  through  the   detector.      (Notice  only  £  =  4,    8,    etc.,   are 
non  zero.  )    t  is  the  "frequency"   relative  to  UU   . 
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Figure  8.      The  first  frequency  block  for  the  LOGjq  amplitude  of  the 
voltage  across  the  detector. 
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Figure  9.      The  first  frequency  block  for  the  relative  phase  between  the 
current  and  voltage  of  the  detector. 
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Figure  10.      The  first  frequency  block  for  the  LOG       amplitude  of  the 
voltage  across  the  junction  without  z    . 
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Figure   11.     The  first  frequency  block  for  the  LOG   Q  power  into  the 
detector. 
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Figure  12.     The  first  frequency  block  for  the  LOG^  power  into  the 
Josephs  on  junction. 
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Figure   13.      The  second  frequency  block  (128  £  l<  256)  for  the  LOGjq 
current  through  the  detector. 
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Figure  14.      The  second  frequency  block  for  the  LOG^  voltage  across  the 
detector. 
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Figure   15.      The  second  frequency  block  for  the  relative  phase. 
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Figure  16.      The  second  frequency  block  for  the  LOG       voltage  across 
junction  without  z   . 
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Figure  17.      The  second  128  lines  for  LOG.q  power  into  detector. 
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Figure  18.     The  second  128  lines  for  LOG       power  out  of  junction. 
(Noise  pulses  are  seen  as  energy  sources  here.  ) 
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(Here  we  have  compared  the  length  of  the  8th  harmonic  to  the  16th  harmonic  which  is  given  as 
full  scale  of  (4.  410).     We  take  the  absolute  value  in  this  effort.     The  minus   sign  is  simply  a  code  to 
tell  us  that  the  power  is  out  instead  of  in.     The  132th  harmonic  has  value  (5.  5  X  (4.  350)/ 6.  2  =  3.  85). 
The  133rd  harmonic  has  value  (3.  5  X  (4.  963)/ 7.  1  =  2.  45).     The  above  values  are  extracted  from  the 
plots  and  the  given  scales  in  the  plots.      To  complete  the  determination  of  the  power  we  need  to  look  at 
the  conversion  table  of  figure  1  9  to  find  that  the  positive  4th  harmonic  is  reduced  by  -2,   the  negative 
8th  is  reduced  by  -6,   the  positive  133rd  is  reduced  by  -24,    and  the  negative  132  is  reduced  by  -16. 
Thus  the  power  is  given  as 


Harmonic 

Power 

4 

103.918 

8 

-lo-2-66 

132 

-10-12'15 

133 

lo"21'55 

To  scale  these  numbers  to  watts  just  multiply  them  by  the  4.  8830  X  10        .     Notice  that  the  noise 
signal  in  the  133  harmonic  is  approximately  9  decades  less  than  the  corresponding  non- noise  signal 
of  the  132  harmonic. 

It  should  be  obvious  that  if  one  wants  extensive  detailed  tables  of  the  various  spectra,    it  is 
better  to  put  in  appropriate  print  statements  into  the  programs.     Our  purpose  has  been  mainly  qualita- 
tive visual  information  of  the  power  spectra.     Figure  19  has  a  table  of  power  and  other  information  at 
selected  harmonics.     These  numbers  are  in  the  dimensionless  form  used  in  the  calculation.     To  con- 
vert to  watts,    volts,    amperes,    etc.  ,    it  is  necessary  to  use  the  constants  given  in  figure  5. 

We  finish  this  section  by  describing  the  remaining  details  shown  in  figure  19.     The  first  row 
of  data  gives  the  total  steady  state  power  across  the  junction  without  the  Z     part  (POWJ),   the  total 
steady  state  power  across  the  full  junction  (POW1),    the  power  across  the  detector  (POW2),    the  power 
put  in  by  the  applied  voltages   (POWIN),    and  finally  the  effective  impedance  KO.     The  next  block  of  six 
different  row  pairs  give  the  harmonic  number  to  be  applied  to  the  row  pair  information  (KEY),   the 
applied  voltages   (B),    the  net  initial  phase  at  s  =  0  (PHC),    the  power  in  (POWIN),   the  power  across 
the  junction  (POW1),    and  the  power  across  the  detector  (POW2).     The  second  row  of  a  given  row  pair 
has  the  dimensionless  values  for  the  complex  number  impedances  of  that  frequency.     The  frequencies 
chosen  to  be  printed  out  are  d.  c.  ,    the  direct  beat  frequency,   the  lower  applied  r.  f.    frequency,   the 
upper  applied  r.  f.   frequency,   the  highest  frequency  used  as  a  consequence  of  the  number  of  data 
points  computed,    and  finally  the  rotation  frequency. 

The  remaining  undiscussed  information  has  the  net  d.  c.  voltage  across,  and  the  current  through  the 
junction  without  the  Z     term,   the  data  summarizing  the  error  controls  used  in  the  program  (EP1, 
EP2,    EP3)  and  finally  the  last  computed  P(a)  set     EP4(I6)    .     Please  look  at  the  program  to  see  what 
these  errors  limits  mean  in  detail. 

6.     THE  SUMMARY 

In  this  report  we  are  documenting  the  computer  program,   the  mathematical  techniques  and 
approximations  used,    and  the  procedures  for  use  and  interpretation  of  the  computer  program.     The 
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Figure  19.     Scaling  information  for  a  particular  steady  state  run. 
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results  using  this  program  are  to  be  shown  in  other  reports. 


Here  the  discussion  is  directed  toward  the  procedures  used  to  check  the  correctness  of  the 

4 


computer  program.     The  strict  d.  c.    cases  have  been  selectively  compared  with  the  Stewart     voltage- 


current  model.      The  r.  f.    applied  voltage  case  has  been  selectively  compared  with  the  Fack,    et  al, 
model  and  their  reported  results.     Both  cases  agree  within  the  accuracy  of  the  techniques  and  the 
graphs  of  the  data  reported  in  those  papers.     The  final  check  has  been  a  comparison  with  the  voltage 
driven  model.     Again  the  results  are  in  agreement.     The  integral  equation  part  has  not  been  checked 
against  any  other  model  or  results  since  thereare  no  such  published  results.     In  terms  of  the  computer 
program,    the  program  has  been  checked  to  be  sure  that  the  D  is  correctly  generated.     The  only 
remaining  part  for  possible  programming  errors  is  in  the  integral  equation  section  in  subroutine  PART. 
The  behavior  and  visual  check  shows  no  obvious  errors.     A  detailed  quantitative  test  is  not  known. 

The  most  important  numerical  conclusion  was  to  establish  that  it  is  very  necessary  to  fix  the 
cp    (o)  and  to  allow  the  V    (o)  to  be  adjusted  when  one  wants  to  pin  the  rotation  frequency.     If  one  wants 
to  use  only  time  domain  information,    then  one  should  use  the  much  simpler  programs  such  as  used  by 
Fack,    Stewart,    and  is  used  in  the  appendix  and  fix  b     instead  of  the  rotation  frequency.     If  one  is 
interested  in  high  order  harmonic  mixing,    and  detailed  harmonic  spectra,   then  this  program  is  best. 
If  noise  effects  on  the  system  of  equations  are  of  interest,   then  it  is  necessary  to  use  the  time  sequence 
procedures  and  to  allow  selected  Fourier  transforms  on  the  various  time  series.     That  analysis  will 
be  discussed  in  another  report. 

Because  the  analysis  has  been  done  for  discrete  uniformly  spaced  time  points,    we  have  discrete 
induced  frequencies  that  are  rationally  related  to  each  other.     The  consequences  of  this  situation  are 
additional  steps  induced  in  the  d.  c.    current-voltage  plots.     The  excursion  of  those  steps  is  reduced 
by  increasing  the  number  of  time  points  used.      To  get  an  idea  on  the  range  of  any  d.  c.    step  for  a  given 
rotation  frequency,    it  is  necessary  to  make  at  least  four  calculations.      First  we  chose  a  particular 
PHA(l)  =  cp     and  compute  the  B     for  that  case.      If  the  program  converges,    then  we  have  one  point  for 
that  case.      To  get  a  convergence  case,    we  use  the  transient  program  with  the  chosen  B     and  deduce  an 
adequate  PHA(l).      The  steady  state  program  just  improves  the  accuracy  of  the  results.     Second,    we 
compute  a  B       with  PHA(l)  now  cp     ±  —  .     At  least  one  of  these  will  converge.     If  both  converge,   then 
(B     +  B    )/ 2  is  the  center  and   |(B      -  B    )|  is  the  range  of  the  step.      For  example,    if  the  B    (B    )  does  not 

T  -  T  -  T  - 

converge  because  the  phase  location  denies  steady  state  solution,    then  it  is  necessary  to  compute 
cp     +  — ,    (cp    -    — )  case  to  get  the  B   _i  ,    (B   j_)  value.     We  assume  that  the  range  of  a  step  is  defined  by  the 
following  simple  formulae--  B  =  X  +  R  sin  9,    where  X  is  the  center  of  the  step,    R  is  the  range  and  9  is 
the  phase  relation  relative  to  the  center  of  step.      Normally  9  is  not  equal  to  cp     and  normally  the  formu- 
lae is  not  strictly  accurate.     It  is  only  true  if  the  dominate  processes  creating  a  step  are  due  to  a 
single  frequencing  mixing.     For  estimating  the  range's  order  of  magnitude,   we  can  accept  this 
approximate  formulae  to  get  the  range  and  center  of  the  step  for  the  case  considered,    namely  B    (B    ),  B 
respectively  are  given  by 

B     =  X  +  R  sin  9 

B_i  =  X  +  R  sin  (9  =F  tt/ 4) 

B     =  X  +  R  sin  (9  =F  tt/2) 
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X=    (B,  +   B^  -"JT  B_i)/(2  -yfT) 

1  +  +2 

R  =  [(B]   -  X)2+   (B^  -  X)2]' 

Bl   "X 

— =  sin  9  or  6  is  gotten  to  give  center.      For  more 

accurate  analysis,    it  is  necessary  to  do  empirical  work. 

7.     APPENDIX 

In  this  section  we  discuss  briefly  a  simpler  computer  program  called  a  transient  program  that 
can  be  used  to  get  improved  estimates  for  the  effects  of  r.  f.    and  d.  c.    applied  voltages  to  the 
Josephson  equations  without  the  operator  D.      The  equations  used  are  closely  related  to  those  given  in 
references   1  and  4.      We  use  the  trapezoidal  spline  method  with  variable  step  size  so  that  the  accuracy 
per    step   is   kept   uniform   throughout.        Figure  20  gives  the  Fortran  listing  for  this  program.      This 
program  simply  generates  a  time  series  of  selected  points.      This  data  can  then  be  used  to  estimate  the 
average  rotation  frequency,    N   ,    and  to  fix  a  PHA(l),    VA(1),    and  a  U. 

7.  1     The  Transient  Equation 

Here  we  show  the  equations  used  in  the  program  and  describe  how  they  relate  to  equation  11. 
The  transient  equations  are: 

*£=  v 

ds 

dV  —      —  — 

—  =    -v^V  -  a  sincp+  B+  BT    cos   (NT  s  +  6T  )+  BTT  cos   (NTTs  +  9TT)  (Al) 

QS  1  Li  Lt  L  U  U  U 

Here  cp  =  cp    +  cp    +  cp     , 

a         b         c 

where  dcp  dV 

HT=  V  ST  =  "vTVb+  *l  cos  (nls  +  9l)+  *u  cos  (V  +  V  ' 

and  cp     is  the  initial  value  of  cp     at  s  =   0;  thus  here  cp     =  0  at  s  =  0.     We  note  that  B  -  b   /K_,    Bt    =  bT/K., 
c  a  a  O       0         Li         L,      0 

and  B      =  b../K   .      The  best  method  to  fix  a,    v_,,    and  K     is  to  use  the  main  program  first  for  a  given 

C,    R,    Z.,    a   ,    J,    etc.  ,    with  a  d.  c.    case.     We  then  use  this  program  called  the  transient  calculation  for 

the  various  B    ,    B    ,    9     ,    9      cases.     Please  note  that  the  scale  used  for  b   ,    h    necessarily  gives 

Li  U  Li  U  Li  U 

different  phases  and  amplitudes  in  the  main  program  if  the  K(N    ),    K(N    )  are  different  from  K   . 

Li  U  \J 

When  the  programs  are  used,   this  difference  needs  to  be  considered.     Before  we  look  at  the  input  and 
output  of  this  program,    we  briefly  consider  the  criteria  for  the  step  size.     It  simply  restricts  the 
absolute  change  in  cp  to  be  less  than  STEP/  (3V1-V2)  and  the  absolute  change  in  V  to  be  less  than 

STEP'W  [vT(3Ql-Q2)].     The  STEP  is  the  control  accuracy  for  this  allowed  change  per  step;  the  VI 

dV 
and  V2  are  the  last  two  v1  s  used;  the  Ql  and  Q2  are  the  last  two  q's  used.     The  q  is  - —  . 

as 
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7.  2    The  Use  of  the  Transient  Program 


This  program  is  organized  under  the  assumption  that  a  given  computer  run  will  fix  the  r.  f.    and 
selected  control  parameters  and  will  sweep  the  applied  d.  c.    current  B  over  many  values.     Thus  the 
first  read  statement  has  a  parameter  which  fixed  the  range  of  s.     The  initial  V,    and  PH  has  no  direct 
correspondence  to  the  expected  output  VA  and  PHA,    hence  they  are  used  to  just  start  the  calculation. 
The  IPRINT  defines  how  many  steps  are  skipped  before  the  usual  print  mode  takes  place.      Figure  21  shows 
the  flow  of  input  data.     Figure  22  gives  an  example  of  printed  output  data.     The  first  three  rows  of 
printed  output  simply  gives  the  summary  of  the  input  data  use  in  this  particular  computer  run.     The 
list  of  names  S,    VA,    etc.  ,    are  the  labels  for  the  six  columns  of  numbers  that  follow  the  Warning 
statement.     The  S  is  the  time  variable.     The  VA  is  the  cp&  with  the  average  steady  state  slope  of  the 
ramp  term  not  removed.     PH  is  the  total  phase  cp.     PHB  is  the  r.  f.    phase  including  the  initial  PHC. 
VR  is  the  integral  of  V;  it  is  used  to  estimate  the  slope  of  the  ramp.     PHA  is  the  cpa  phase  with  the 
ramp  term  still  present.      To  illustrate    how  to  get  the  necessary  initial  data  for  the  main  program 
from  this  transient  program,   we  run  through  the  sample  output  to  get  the  proper  numbers.     Please 
note  that  the  printed  output  far  exceeds  the  necessary  information  for  these  variables.     It  may  be 
possible  to  suppress  all  printing  except  the  few  numbers  of  interest.     The  current  calculation  assumes 
that  the  transients  are  unimportant  after  S  exceeds  2tt,    hence  we  let  the  computation  only  run  about  one 
more  2tt  interval.     This  last  interval  is  used  to  estimate  the  ramp  and  the  needed  initial  PHA(l),    VA(1), 
and  U  in  the  main  program.     Note  that  Pl6,    Vl6,    and  BI  are  respectively  used  as  the  input  variables 
for  this  initial  data.     In  our  sample  output  data  we  have  the  following  needed  data: 

s  VA  PHB  VR  PHA 

6.28245  8.665  -.756  20.74  21.486 

12.  57033  8.933  -.645  46.03  46.66 

B  =   102,    and  v      =  32.     We  use  subscript  1  and  2  to  label  each  data  now.     We  chose  these  points 
because  we  want  PHB  to  be  nearly  the  same  in  both  this  program  and  in  the  main  program.     If  you 
change  to  other  THL  and  THU  then  you  can  consider  other  data  points.     Keep  in  mind  that  it  is  not 
necessary  to  have  accurate  values.     The  purpose  is  to  get  a  guess  that  is  at  least  accurate  to  10  to  20% 
so  that  the  circle  of  convergence  for  the  quasilinear  procedure  will  operate  effectively  in  the  main 
program  and  so  that  some  information  is  available  to  estimate  potential  effects  due  to  spiking.     We 
can  then  choose  the  necessary  accuracy  in  the  main  calculation.     The  ramp  slope  is  estimated  as 


N. 


J    =   (VR2  -  VRj)/(s2  -  Bj)  +  4. 
The  initial  data  is  thus  estimated  as 


0 


BI  =  B  -  VT*N    =  -26 


VI 6  =   (VA    +  VA2)/2.    -  Nj  =  4.  8 
PI  6  =      (PHA     -  N  *2tt)+  (PHA2  -  Nj*4n)    /  2. 
=  (PHA    +  PHA  /2.   -  N  *3tt  =  -3.626  . 

1  Ct  I 
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PROGRAM  JOSPA 
C      TO  GET  INTIAL  WAVE  SHAPE  FOR  STARTING  GUESSES  IN  JOSPHC 
C      PHDOT=V   $  VDOT=B-V*ANUT-A*SINF(PH>+BL*COSF<NL*S+THL>+ 
C      BU*COSF(NU*S+THU) 

COMMON  /l/  B.ANUT.BL.THL.BU.THU.NL.NU.A.STEP.STEP1 

READ  1.RANGE»V.PH,BL.BU,THL.THU»STEP.ANUT,A.IPRINT,NL.NU 

70  READ  l.B 

1  F0RMAT(5E12.3/5E12. 3/3110) 
IF(EOF»60)71»72 

71  CALL  EXIT 

72  PRINT  2.B.V.PH.BL.BU.THU.THL.STEP.ANUT.A.RANGE.IPRINT  »NL»NU 

2  F0RMAT(*  B.V.PH.BL.BU/THU.THL.STEP.ANUT.A.RANGE/IPRINT  .NL,Nu  * 
1  /5E12.3/6E12. 3/3110/*  S.VA.PH.PHB.  VR*PHA  *  //) 

PRINT  99 
99     FORMAT <*  WARNING  THE  VA  AND  PHA  HAVE  NOT  HAD  THE  RAMP  TERMS  REMOVE 
ID  *) 

J=0  $STEP1=STEP*A/ANUT$IF(ABSF(V)-1.)26.26.25 

26  H=STEP/100.  $GO  TO  270 
25     H=STEP/( 100.*ABSF(V) ) 
270    PHC=PH         $  S=0. 

Il=-2$RANGl=RANGE-7. 

27  11=11+1   $  IF<S-RANGE)4,4»70 

40  J=J-1         $IF(J)41»41»3 

41  IF(S-RANG1>3. 417.417 
417  IF(S)410»410,412 
410  IF<NL*NU>70»70»415 

415         Xl=-BL/(NL*SQRTF(NL*NL+ANUT*ANUT) ) 
X2=-BU/(NU*SQRTF<NU*NU+ANUT*ANUT) ) 
APH1=ATANF( ANUT/ND       +THL      $X3=-X1*NL 
APH2=ATANF(ANUT/NU)       +THU      $    X4=-X2*NU 

412  S8=NL*S+APH1      $    S9=NU*S    +APH2 
PHB=X1 *COSF ( S8  >  +X2*C0SF ( S9  > 
VB=X3*SINF(S8)+X4*SINF(S9) 
VA=V-VB 

PHA=PH-PHB  -PHC 

413  PRINT  42.S.VA.PH.PHB.VR  »PHA  $J=IPRINT 

3  CONTINUE   S  GO  TO  27 

42  F0RMAT(7E15«8) 

4  IF( Il>14»15»16 

14  V2=V  $  PH2=PH  $S=S2=0.  $  CALL  START < PH. S»QX )  $OT=QX-ANUT*V  $60  TO 
1   410 

15  PHY=PH2+H*V2  SCALL  COM ( PH.V.Q .S .PH2 .V2 .QX .0. .H.PHY) 
VR=V*H 

QT=Q-ANUT*V  £  PH1=PH$  V1=V  $Q2=QX  $Q1=Q  $S1=S  $GO  TO  40 

16  CALL  PR0J(PHY,H1.H»V2»V1»Q2»Q1»PH1) 

CALL  C0M(PH.V.Q,S,PH1»V1.Q1.S1.H1»PHY)  $QT=Q-ANUT*V 

DL1=PH-PH1  $DL2=PH1-PH2 

Q2=Q1   $V2=V1  $PH2=PH1  $S2=S1 

VR=VR+V*H1 

Q1=Q  $  V1=V  $  PH1=PH  $S1=S  $H=H1  $IF (DL1*DL2 ) 412 .412 .40 

END 


SUBROUTINE  COM(PH,V»Q.S.PHX . VX »QX .SX.HX.PHY) 

COMMON  /l/  B.ANUT.BL.THL.BU.THU.NL.NU.A.STEP.STEPl 

DATA(ERO=l.E-5) 

H=HX 

H3=H*ANUT/2.  $H4=1.-H3  $H5=1.+H3  $S=SX+HX 

CALL  QR(O.QY.QZ.PHY)  $  H6= { H*H*.25 ) /H5  $H9=H4/H5  $H10=H/ ( 2.*H5 ) 

H7=1./<1.-H6*QZ)  $H8=VX*H/H5+PHX+H6*QX  $QZ1=QF(S)  $QY=QY+QZ1 

IV=5 

4  IV=IV-1  $  IF( IV)5.6»6 

6  DEL=(H8-PHY+H6*QY)*H7 

PHY=PHY+DEL    SCALL   QR ( 1 »QY.O. .PHY )    SQY=QY+QZ1 
IF(ABSF(DEL>-ERO)8»4»4 

5  PRINT    7.S    SCALL    EXIT 

7  FORMAT (*    NO    CONVERGENCE    *.E12.3) 

8  PH=PHY    $Q=QY    $V=(H9*VX+H10*(QY+QX)     )  $RETURN    $END 


Figure  20.      The  transient  Fortran  program  listing. 
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SUBROUTINE  PROJ (PHY. HI ,H »V2 » VI ,Q2 *Ql »PH1 ) 

COMMON  /l/  B»ANUT,BL»THL»BU»THU.NL.NU»A*STEP.STEPl 

H1A=STEP/ABSF(3.*V1-V2) 

H1B=STEP1/ABSF(3.*Q1-Q2> 

H1=MIN1F(H1A,H1B) 

PHY=PH1+H1*(VH-(V1-V2)*H1/(2.*H) )       SRETURN  SEND 


SUBROUTINE  START ( PHX »SX »QX ) 

COMMON  /l/  B.ANUT,BL»THL»BU»THU»NL»NU»A»STEP»STEP1 

CALL  QR( 1»QX»0. .PHX)  $QX=QX+QF ( SX >  SRETURN  SEND 


FUNCTION  QF(Z) 

COMMON  /l/  B.ANUT.BL.THL.BU.THU.NL.NU.A.STEP.STEP1 

QF=B+BL*COSF(NL*Z+THL)+BU*COSF(NU*Z+THU)  SRETURN  SEND 


SUBROUTINE  QR( I .QA.QB.X) 

COMMON  /l/  B»ANUT.BL»THL»BU»THU»NL»NU»A.STEP»STEPl 

IF(  Dl»l»3 
1  QB=-A*COSF(X) 

3  QA=-A*SINF(X)     SRETURNSEND 


Figure  20.    (Continued) 
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PROGRAM  JOSPA 

READ  1 ,  RANGE,  ETC 

70  READ  1 ,  B 
IF(EOF)  EXIT 
IF(S-RANGE)4,4,70 


COMPLETE 
PRINT 


SUBROUTINES    CALLED 

COM 

PROJ 

START 

QF 

QR 


Figure  21.      The  flow  chart  for  Josephson  transient  program. 
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Note  that  only  integer  ft     can  be  used  in  the  main  program  and  that  the  choice  of  B  in  the  transient 
program  does  not  necessarily  correspond  to  such  integers. 

8.     ACKNOWLEDGMENTS 

We  wish  to  thank  Stephen  Jarvis  and  Clark  Hamilton  for  useful  discussions. 

9.     REFERENCES 

1.  W.    C.   Stewart,  App.   Phys.   Letters  1_2,    8(1968),   pp.    277-280. 

2.  Theory  and  Application  of  Spline  Functions,    Ed.    T.   Greville,   Academic  Press  (1969).     An 
Introduction  to  the  Application  of  Spline  Functions  to  Initial  Value  Problems,   Frank  R.   Loscalzo, 
pp.   37-64. 

3.  A.   Miele,   "Modified  Quasilinearization  Method  for  Solving  Nonlinear,   Two  Point  Boundary  Value 
Problems,"  (Rice  University,    1970),    Clearinghouse  Reports  N71-1779. 

4.  H.   Fack,   V.   Kose,   andH. -J.   Schrader,   "Calculation  of  Characteristic  Curves  for  Point- Contact- 
Josephson- Elements  Subjected  to  High  Frequency  Radiation,"  Messtechnik  (February  1971)  pp.    31-38. 


55 


FORM   NBS-114A    (1-71) 


U.S.   DEPT.   OF  COMM. 

BIBLIOGRAPHIC  DATA 
SHEET 


1.  PUBLICATION  OR  REPORT  NO. 

NBS  TN-627 


2.  Gov't  Accession 
No. 


3.  Recipient's  Accession  No. 


4.  TITLE  AND  SUBTITLE 

Computation  of  Spectral  Data  for  a  Josephson 
Junction  Circuit 


5.    Publication  Date 

November  1972 


6.  Performing  Organization  Coc 


7.  AUTHOR(S) 

E.    G.    Johnson,    Jr.    and  D.    G.    McDonald 


8.  Performing  Organization 


10.  Project/Task/Work  Unit  N, 

2710164 


9.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

NATIONAL  BUREAU  OF  STANDARDS,     Boulder   Labs 
DEPARTMENT  OF  COMMERCE 

Boulder,    CO    80302 


11.  Contract/Grant  No. 


12.  Sponsoring  Organization  Name  and  Address 


Same  as  9  above. 


13.  Type  of  Report  &  Period 
Covered 

Final  (4/72-10/72) 


14.  Sponsoring  Agency  Code 


15.  SUPPLEMENTARY  NOTES 


16.  ABSTRACT  (A  200-word  or  less  factual  summary  of  most  significant  information.    If  document  includes  a  significant 
bibliography  or  literature  survey,  mention  it  here.) 

A  computer  program  has  been  developed  to  study  power  flow  between  differem 
frequency  channels  in  a  Josephson  junction  circuit.  This  paper  discusses  the 
mathematical  assumptions  used  to  get  such  results.  They  are  the  trapezoidal 
approximation  from  spline  theory  and  the  use  of  a  finite  range  of  frequencies 
characterize  the  frequency  spectrum.  This  paper  describes  the  program  and 
provides  the  FORTRAN  listing,  flow  charts,  and  discusses  how  to  use  the 
program.     A  discussion  of  possible  sources  of  errors  is  also  included. 
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